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ABSTRACT 


This  report  presents  a  2-D  depth- averaged  numerical  model  simulating  water  flow  and 
reactive  contaminant  and  sediment  transport  in  bay/estuary  systems.  This  model  includes  two 
computational  modules:  flow  and  transport.  We  use  a  splitting  strategy  to  solve  hydrodynamic 
equations  where  the  method  of  characteristics  (MOC)  is  employed  first  to  compute  circulation 
without  taking  into  account  eddy  fluxes,  and  then  the  Galerkin  finite  element  method  is  applied  to 
deal  with  eddy  fluxes.  The  computed  flow  results  are  used  to  compute  the  migration  of  sediments 
and  contaminants  in  the  transport  module.  Sediments  can  be  either  suspended  in  the  water  column 
or  deposited  on  bed  rock,  and  whose  distributions  are  determined  subject  to  advection,  dispersion, 
and  deposition/erosion  processes.  Contaminants,  either  organic  or  inorganic,  may  exist  in  the 
dissolved  phase  in  the  water  column  or  in  the  interstitial  water  of  the  bed  sediments.  They  may 
appear  also  as  particulate  chemicals  on  sediments  through  adsorption  reactions.  All  chemical 
reactions,  including  aqueous  complexation,  adsorption/desorption,  and  volatilization,  are  assumed 
elementary  and  kinetic  with  reaction  rates  calculated  based  on  the  collision  theory  where  forward 
and  backward  rate  constants  are  determined  experimentally.  We  employ  a  predictor-corrector 
strategy  in  the  transport  module  where  transport  equations  are  computed  in  the  predictor  step  with 
reaction  rates  estimated  at  the  previous  time,  and  rate  correction  is  achieved  node  by  node  in  the 
corrector  step  to  obtain  the  concentration  distributions  at  the  present  time.  The  Lagrangian-Eulerian 
approach  is  applied  to  solve  the  transport  equations  in  the  predictor  step  with  a  particle  tracking 
technique  used  to  handle  advection.  The  Picard  method  and  the  Newton-Raphson  method  are  used 
to  solve  the  rate  correction  equations  for  sediment  distribution  and  reactive  chemistry,  respectively, 
in  the  corrector  step.  Four  examples  are  used  to  demonstrate  the  capability  of  this  model. 
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Chapter  1.  INTRODUCTION 

For  at  least  two  decades,  environmental  issues  have  become  major  public  concerns 
worldwide  and  have  drawn  scientists  and  engineers  from  a  variety  of  disciplines  to  devote 
themselves  to  environmental  impact  analysis/assessment  (EIA)  studies.  Along  with  the  rapid 
development  of  computer  technology  and  the  advancement  of  knowledge  in  numerical  methods, 
numerical  simulation  techniques  have  become  extremely  useful  tools  in  today’s  EIA  studies.  Many 
numerical  models  have  been  presented  to  simulate  circulation  of  fluid  flow  and/or  mass  transport 
of  contaminants  and  sediments  in  order  to  deal  with  problems  relating  to  both  water  quantity  and 
quality  in  bay /estuary  areas  [Falconer  and  Lin,  1997;  Paulsen  and  Owen,  1996;  Blumberg  et  al., 
1996;  Falconer,  1990;  Lai,  1977], 

Many  numerical  strategies  have  been  presented  to  handle  the  non-linear  hydrodynamic  flow 
equations  that  are  used  to  describe  the  circulation  of  estuarine  waters.  The  numerical  strategies 
included  total  variation  diminishing  method  [Tseng  and  Liang,  1997],  two-step  explicit  method 
[Park  et  al.,  1995],  semi-implicit  method  [Cecchi  et  al.,  1998],  implicit  method  [Fennema  and 
Chaudhry,  1989],  high-order  time  integration  schemes  [Ozkan-Haller  and  Kirby,  1997], 
upwind/upstream  numerical  schemes  [Shi  and  Toro,  1996;  Anastasiou  and  Chan,  1997],  filtered 
solution  method  [Laible  and  Lillys,  1997],  least  square  method  [Muccino  et  al,  1997;  Tseng  and 
Liang,  1997],  the  method  ofcharacteristics  [Lai,  1977;  Hirsch  et  al.,  1987],  and  hybrid  Lagrangian- 
Eulerican  (LE)  approach  [Hansen  and  Arneborg,  1997;  Petera  and  Nassehi,  1996].  On  the  other 
hand,  a  variety  of  numerical  strategies,  such  as  splitting  method  [Falconer  and  Lin,  1997;  Sommeijer 
andKok,  1996],  upwind  fomulation  [Falconer  and  Lin,  1997;  Ors,  1997;  Sommeijer  and  Kok,  1996]; 
LE  approach  [Spitaleri  and  Corinaldesi,  1997;  Wood  and  Baptista,  1992],  and  decoupling  scheme 
[Park  and  Kuo,  1996],  have  been  introduced  to  handle  contaminant  and  sediment  transport  equations. 

In  the  past  decade,  the  LE  approach  has  attracted  many  numerical  modelers’  attention 
because  it  can  avoid  many  types  of  numerical  errors  when  advection/convection  dominates  [Yeh, 
1990].  Casulli  and  Cheng  have  shown  that  the  approach  is  unconditionally  stable  when  it  was  used 
to  solve  1-D  shallow  water  equations  [Casulli  and  Cheng,  1990].  Petera  and  Nassehi  employed  a 
particle  tracking  technique  in  the  Lagrangian  step  and  have  demonstrated  their  approach  would 
provide  very  stable  results  for  the  convection  dominated  very  shallow  depth  computations  in 
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estuaries  [Petera  and  Nassehi,  1996].  Yeh  and  his  colleagues  have  proposed  adaptive  local 
refinement  numerical  algorithms  to  accurately  solve  advection-diffusion  transport  equations  on  the 
LE  framework  [Cheng  et  al.,  1996a;  Yeh  et  al.,  1992].  In  our  2-D  depth-averaged  model  presented 
here,  we  also  used  the  LE  approach  to  handle  the  advection  and  dispersion  terms  in  both  flow  and 
transport  equations. 

The  method  of  characteristics  (MOC)  might  be  considered  the  most  adequate  way  to  solve 
wave  or  hyperbolic-type  equations  because  it  approaches  problems  not  only  mathematically  but 
physically  also.  Lai  has  developed  comprehensive  MOC  models  to  solve  1-D  shallow  water 
equations  [Lai,  1987]  and  has  incorporated  the  MOC  with  the  finite  difference  discretization  to  solve 
2-D  flow  equations  [Lai,  1977].  In  this  work,  we  apply  the  MOC  to  solve  flow  equations  on  the 
platform  of  finite  element  discretization.  Diagonalization  [Abbott,  1966]  can  be  achieved  easily 
when  the  MOC  is  used  to  solve  1-D  equations.  Its  application  to  2-D  cases,  however,  is  not 
straightforward  [Hirsch  et  al.,  1987].  The  three  characteristics  associated  with  2-D  shallow  water 
equations  can  be  chosen  in  any  direction  [Lai,  1977].  Many  iterations  might  be  needed  to  reach 
convergence  in  the  numerical  simulation  if  the  characteristic  directions  are  not  appropriately 
selected.  To  handle  this,  we  perform  local  diagonalization  for  the  2-D  characteristic  equations.  We 
also  apply  the  Lagrangian  approach  to  solve  the  characteristic  equations,  rather  than  the  original 
governing  equations,  where  a  particle  tracking  technique  [Cheng  et  al.,  1996b]  was  employed  to 
transport  characteristic  variables  along  characteristics. 

Though  many  models  have  included  updated  numerical  approaches  to  improve  their 
computational  performance,  no  existing  water  quality  model,  to  our  knowledge,  has  used  a  generic 
and  mechanistic  approach  to  simulate  the  reaction-migration  behavior  of  non-conservative  chemicals 
in  bay /estuary  areas .  All  existing  models  are  either  empirically-based  or  simulate  systems  containing 
specific  reactions,  and  thus  their  applications  are  limited  to  specific  systems.  Unlike  an  empirically- 
based  model  whose  model  parameters  might  be  system-  or  even  environment-dependent,  a  generic 
and  mechanistic  model  considers  interactions  among  chemicals  based  on  reaction  mechanisms  which 
are  universal  and  can  be  applied  to  all  systems.  Based  on  our  experience  of  modeling  reactive 
contaminant  transport  in  the  subsurface  [Cheng  and  Yeh,  1998;  Yeh  et  al.,  1998;  Yeh  and  Tripathi, 
1990],  we  use  a  generic,  mechanistic  approach  to  simulate  reactive  contaminant  transport  in 
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bay /estuary  areas,  desiring  to  enhance  reactive  contaminant  transport  modeling  in  this  field.  In  our 
transport  module,  we  take  into  account  suspended  sediments  and  bed  sediments  whose  distributions 
are  determined  through  hydrological  transport  and  erosion/deposition  processes.  Any  contaminant 
considered  in  our  model  is  treated  as  a  chemical  species  which  can  be  either  organic  or  inorganic. 
It  may  appear  in  the  water  column  as  a  dissolved  chemical,  be  adsorbed  onto  suspended  or  bed 
sediments  to  become  particulate  chemicals,  or  exist  in  the  interstitial  water  of  the  bed  sediments  also 
in  dissolved  form.  All  reactions  among  chemicals  are  assumed  elementary  kinetic,  and  reaction  rates 
are  described  based  on  the  collision  theory  [Stumm  and  Morgan,  1981]  where  forward  and  backward 
rate  constants  can  be  measured  experimentally.  Figure  1.1  shows  a  schematic  plot  for  the  chemical 
reactions  taken  into  account  in  our  model. 


Water  Surface 


Suspended  Sediment 


Water 


Bed  Rock/Soil 


(!)(!’)  Aqueous  complexation  reations  in  the  water  column  and  in  the  interstitial  water 

(2)  Adsorption! desorption  between  the  water  column  and  suspended  sediments 

(3)  Adsorption/desorption  reactions  between  the  water  column  and  bed  sediments 

(4)  Adsorption/desorption  reactions  between  the  interstitial  water  and  bed  sediments 

(5)  Volatilization  reactions  between  the  water  column  and  atmosphere 


Figure  1.1.  Schematic  plot  of  chemical  reactions  taken  into  account  in  the  transport  module. 


In  developing  this  model,  we  have  assumed  (1)  flow  fields  are  not  influenced  by  sediment 
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or  contaminant  distributions  and  (2)  sediment  distributions  are  not  affected  by  contaminant 
distributions.  The  first  assumption  lets  our  model  compute  flow  and  transport  module  sequentially, 
which  as  a  result  demands  much  less  computer  time  than  that  accounting  for  strong  coupling 
between  flow  and  transport  modules  due  to  density  effect  [Cheng  and  Yeh,  1998].  The  second 
assumption  allows  us  to  calculate  sediment  transport  prior  to  contaminant  transport  in  our  numerical 
approach  as  will  be  stated  later  on.  Theoretically,  strong  coupling  [Cheng  and  Yeh,  1998]  must  be 
used  to  couple  chemical  reactions  with  physical  transport  to  provide  accurate  solutions.  This  may, 
however,  require  many  coupling  iterations  to  reach  convergence,  which  might  cause  unaffordable 
computer  effort  for  long-term  simulation  in  surface  systems  where  the  time-step  size  is  in  the  order 
of  minute  or  second.  To  achieve  efficient  computations,  most  reactive  contaminant  transport  models 
employ  operator  splitting  strategies  to  decouple  solute  transport  from  chemical  reactions  [Salvage, 
1998;  Chilakapathi,  1995].  Rather  than  operator  splitting,  we  use  a  predictor-corrector  scheme  to 
effectively  solve  reactive  contaminant  transport  without  surrendering  much  accuracy.  We  also 
employ  the  LE  approach  to  deal  with  transport  equations  where  backward  particle  tracking  is  applied 
in  the  Lagrangian  step  and  the  Galerkin  finite  element  method  is  implemented  in  the  Eulerian  step. 

To  describe  our  model  in  detail  in  the  following,  we  will  first  describe  the  governing 
equations  and  the  associated  initial  and  boundary  conditions  in  Chapter  2.  Then  we  will  state  the 
numerical  approaches  in  Chapter  3,  and  use  an  example  problems  for  demonstration  in  Chapter  4. 
Finally,  we  will  summarize  in  Section  5. 
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Chapter  2.  MATHEMATICAL  BASIS 

In  this  chapter,  we  are  to  give  governing  equations,  initial  conditions,  and  boundary 
conditions  for  simulating  bay  estuary  circulation  as  well  as  contaminant  and  sediment  transport. 

2.1  Flow  Equations 

The  governing  equations  of  2-D  circulation  flow  in  the  bay/estuary  area  can  be  described  by 
shallow  water  equations  which  are  derived  based  on  the  conservation  law  of  water  mass  and  linear 
momentum  and  through  the  depth-averaging  process  [Wang  and  Connor,  1975].  The  governing 
equations  of  a  dynamic  wave  model,  without  taking  into  account  Coriolis  forces,  can  be  written  as 
follows. 

The  continuity  equation: 

3n  5uh  dvh  c 

5t  dx  dy 

where  T)  =  T|(x,y,t)  is  water  surface  elevation  with  respect  to  mean  sea  level;  h(x,y,t)  is  water  depth 
[L];  u  is  the  x-velocity  [L/T];  v  is  the  y-velocity  [L/T];  S  represents  the  external  source  [L/T].  Here 
we  also  define  sea  bed  elevation  Z0  =  -d(x,y),  where  d(x,y)  is  water  depth  at  mean  sea  level  (Figure 
2.1).  It  is  thus  obvious  that  T)(x,y,t)  =  Z0(x,y)  +  h(x,y,t). 


Figure  2.1.  Schematic  definitions  of  water  stage  T)(x,y,t),  water  depth  h(x,y,t), 

and  sea  bed  elevation  Z0(x,y). 

The  momentum  equations: 
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dx  dx 


gh-^-  + 


dF 


xy 


dy  dx 


where  (us,vs)  is  projected  velocity  of  the  external  source  on  the  horizontal  plane  [L/T];  p  is  fluid 
density  [M/L3] ;  g  is  the  acceleration  due  to  gravity  [L/T2] ;  (Txs,Tys)  is  the  surface  shear  stress  per  unit 
horizontal  area  [M/L/T2];  (Txb,Tyb)  is  the  bottom  shear  stress  per  unit  horizontal  area  [M/L/T2]  that 
can  be  estimated  with  Manning’ s  formula  [Chow,  1973] ;  Fxx  and  Fyx  are  the  fluxes  due  to  turbulence 
and  diffusion  along  the  x-direction,  while  Fxy  and  Fyy  are  along  the  y-direction  [L3/T2] .  Fxx,  Fyy,  Fxy, 
and  Fyx  are  hypothesized  to  be  proportional  to  both  water  depth  and  the  derivatives  of  velocity, 
described  as  follows. 
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XX 


he. 


5u 

dx 
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yy 


dv 

dy 


du  +  d 
dy  d 


where  exx,  e  ,  and  exy  are  eddy  coefficients  [L2/T].  To  perform  transient  simulations,  water  surface 
elevation  and  velocity  must  be  given  as  the  initial  conditions.  From  a  physical  point  of  view,  three 
types  of  boundary  conditions  may  be  taken  into  account,  described  as  follows. 

Water  surface  elevation-specified  boundary  condition: 

This  condition  is  applied  when  the  water  surface  elevation  at  a  boundary  node  can  be 
prescribed  through  the  simulation  period.  Such  a  condition  can  be  expressed  as 

Tl(xb,yb,t)  =  T]b(t) 


where  (xb,  yb)  is  the  coordinate  of  a  boundary  node;  T|b(t)  is  a  prescribed  time-dependent  water 
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surface  elevation  [L] . 

Normal  flux-specified  boundary  condition: 

This  condition  is  employed  when  the  time-dependent  incoming  normal  flux  at  a  boundary 
segment  can  be  prescribed.  It  can  be  expressed  as 

)  =  -n-V(xb,yb,t)h(xb, 

where  n  is  the  outward  unit  vector  normal  to  the  boundary  segment;  q  is  the  flow  flux  vector  [L2/T] ; 
V  =  (u,  v);  qb(t)  is  a  time-dependent  normal  flux  [L2/T]. 

Impermeable  boundary  condition: 

This  condition  is  applied  when  the  boundary  segment  is  a  closed  boundary  through  which 
no  water  flow  is  allowed.  It  can  be  expressed  as 

,»t)  =  -n-V(xb,yb,t)h(x 


2.2.  Contaminant  and  Sediment  Transport  Equations 

The  2-D  depth-averaged  governing  equations  of  transport  in  bay/estuary  areas  can  be  derived 
according  to  the  conservation  law  of  material  mass  and  through  the  averaging  process  over  water 
depth  [Wang  and  Connor,  1975].  In  the  model,  suspended  sediments,  dissolved  chemicals  in  the 
water  column,  and  particulate  chemicals  on  suspended  sediments  are  subject  to  transport,  while  bed 
sediments,  particulate  chemicals  on  bed  sediments,  and  dissolved  chemicals  in  the  interstitial  water 
are  taken  as  immobile.  Their  governing  equations  can  be  written  as  follows. 

Continuity  equations  for  suspended  sediments:  n6  [1,NJ 


isn)  = 


hK-VS„ 


+  M 


where  Ns  is  the  number  of  sediment  sizes;  h  is  water  depth  [L] ;  t  is  time  [T];  V  is  the  del  operator; 
Sn  is  depth- averaged  sediment  concentration  of  the  n-th  fraction  size  in  the  water  column  [M/L3];  K 
is  dispersion  coefficient  tensor  [L:/T] ;  Mns  is  source  of  the  n-th  size  fraction  sediment  [M/L2/T] ;  Rn 
is  erosion  rate  of  the  n-th  size  fraction  sediment  [M/L2/T];  Dn  is  deposition  rate  of  the  n-th  size 
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fraction  sediment  [M/L2/T] .  The  determination  of  Dn  and  Rn  can  be  found  in  sediment  transport  texts 
[Graf,  1984], 

Continuity  equations  for  bed  sediments:  n 6  [1,NJ 


at 


-  D  -R 


n 


where  Mn  is  sediment  mass  per  unit  horizontal  bed  area  of  the  n-th  size  fraction  [M/L2]. 
Continuity  equations  for  dissolved  chemicals  in  the  water  column:  ie  [1,NC] 


k; 


af 


hCjW  +  hkjCpj  -  — ^Cj1 

k. 


»  +  S  (amj-bmj)kmh 


m=  1 

\  / 


N 


J  = 


c;w  + 


N, 


n=l  Pn  j 


where  Nc  is  the  number  of  dissolved  chemicals;  C,w  is  depth-averaged  concentration  of  the  i-th 
dissolved  chemical  in  the  water  column  [M/L3];  M™  is  source  of  the  i-th  dissolved  chemical 
[M/L2/T];  A,w  is  combined  first  order  degradation  rate  constant  of  the  i-th  dissolved  chemical  [1/T]; 
kiab  and  k;af  are  backward  and  forward  volatilization  rate  constants,  respectively,  associated  with  the 
i-th  dissolved  chemical  in  the  atmosphere;  p,  is  the  partial  pressure  in  the  atmosphere  associated  with 
the  i-th  dissolved  chemical  [atm] ;  knisb  and  knist  are  backward  and  forward  adsorption  rate  constants, 
respectively,  associated  with  the  i-th  particulate  chemical  on  suspended  sediment  of  the  n-th  fraction 
size;  Cnib  is  mass  of  particulate  chemical  on  unit  mass  of  bed  sediment  of  the  n-th  fraction  size 
[M/M];  knibb  and  knibl  are  backward  and  forward  adsorption  rate  constants,  respectively,  associated 
with  the  i-th  particulate  chemical  on  bed  sediment  of  the  n-th  fraction  size;  Cnis  is  mass  of  particulate 
chemical  on  unit  mass  of  suspended  sediment  of  the  n-th  fraction  size  [M/M];  Nrx  is  the  number  of 
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aqueous  complexation  reactions;  kmrb  and  kmrl  are  backward  and  forward  rate  constants,  respectively, 
of  the  m-th  aqueous  complexation  reaction;  amj  (bmj)  is  stoichiometric  coefficient  of  the  j-th  dissolved 
chemical  in  the  m-th  aqueous  complexation  reaction  when  this  dissolved  chemical  appears  as  a 
reactant  (product)  chemical;  Cibw  is  the  concentration  of  the  i-th  dissolved  chemical  in  the  interstitial 
water  [M/L3];  0n  is  porosity  of  the  n-th  size  fraction  bed  sediment;  pn  is  bulk  density  of  the  n-th  size 
bed  sediment  [M/L3];  E  is  exchange  coefficient  due  to  diffusion  between  pore  water  and  column 
water  [L3/L2/T] ,  which  can  be  determined  experimentally.  In  Eq.  (2.12),  the  first  term  on  the  left  side 
represents  the  mass  increase  rate  of  a  dissolved  chemical  per  unit  horizontal  area,  while  the  second 
term  on  the  left  is  the  advection  transport  term.  On  the  right  hand  side  of  Eq.  (2.12),  the  first  term 
is  the  dispersion  transport  term,  the  second  is  external  source,  the  third  represents  an  overall  first 
order  decay,  the  fourth  is  from  the  volatilization  reactions,  the  fifth  and  the  sixth  are  from  the 
reactions  of  adsorption/desorption  between  the  water  column  and  suspended  sediments  and  between 
the  water  column  and  bed  sediments,  respectively,  the  seventh  is  from  the  aqueous  complexation 
reactions  in  the  water  column,  the  eighth  is  the  exchange  due  to  diffusion  between  the  water  column 
and  the  interstitial  water,  the  ninth  and  tenth  are  the  physical  exchanges  between  the  water  column 
and  the  interstitial  water  due  to  deposition  and  erosion,  respectively. 

Continuity  equations  for  particulate  chemicals  on  suspended  sediments:  n£  [1,NJ,  ie  [1,NC] 

cn-) 

)l+M:;-x:1hsicn;-k 


where  Mnjcs  is  source  of  the  i-th  particulate  chemical  on  suspended  sediment  of  the  n-th  fraction  size 
[M/L2/T];  Anis  is  combined  first  order  degradation  rate  constant  of  the  i-th  particulate  chemical  on 
suspended  sediment  of  the  n-th  fraction  size  [1/T].  In  Eq.  (2.13),  the  first  term  on  the  left  side 
represents  the  rate  of  mass  increase  of  a  suspended  particulate  chemical  per  unit  horizontal  area, 
while  the  second  term  on  the  left  is  the  advection  transport  term.  On  the  right  hand  side,  the  first 
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term  is  the  dispersion  transport  term,  the  second  is  external  source,  the  third  represents  an  overall 
first  order  decay,  the  fourth  is  from  the  reactions  of  adsorption/desorption  between  the  water  column 
and  suspended  sediments,  the  fifth  is  due  to  erosion,  and  the  last  is  due  to  deposition. 

Continuity  equations  for  particulate  chemicals  on  bed  sediments:  nE  [1,NJ,  iE  [1,NC] 


R„c„“) 


k2f»  r  bw 

n,MnCi 


where  Anib  is  combined  first  order  degradation  rate  constant  of  the  i-th  particulate  chemical  on  bed 
sediment  of  the  n-th  fraction  size  [1/T].  In  the  above  equation,  the  term  on  the  left  side  represents 
the  rate  of  mass  increase  of  a  bed  particulate  chemical  per  unit  horizontal  area.  The  first  and  second 
terms  on  the  left  of  the  equation  are  due  to  deposition  and  erosion,  respectively,  the  third  represents 
an  overall  first  order  decay,  the  fourth  is  from  the  reactions  of  adsorption/desorption  between  the 
water  column  and  bed  sediments,  and  the  fifth  and  sixth  are  from  the  reactions  of  adsorption/ 
desorption  between  the  interstitial  water  and  bed  sediments. 

Continuity  equations  for  dissolved  chemicals  in  the  interstitial  water:  nE  [1,NJ,  iE  [1,NC] 


where  kni2b  and  kni21  are  backward  and  forward  adsorption  rate  constants,  respectively,  associated  with 
the  i-th  particulate  chemical  in  the  interstitial  water  of  bed  sediment  of  the  n-th  fraction  size.  In  Eq. 
(2.15),  the  term  on  the  left  side  represents  the  mass  increase  rate  of  a  dissolved  chemical  in  the 
interstitial  water  per  unit  horizontal  area.  The  first  term  on  the  left  is  due  to  diffusion  between  the 
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water  column  and  the  interstitial  water,  the  second  is  from  deposition,  the  third  is  from  the  reactions 
of  adsorption/desorption  between  the  water  column  and  bed  sediments,  and  the  last  is  from  the 
aqueous  complexation  reactions  in  the  interstitial  water. 

In  performing  Eqs.  (2.12)  through  (2.15),  the  following  reactions  are  considered. 

Aqueous  complexation  reactions  in  the  water  column  and  in  the  interstitial  water:  m6  [l,Nrx] 

.  T* 

1 

i=l  i=  1 

k*fl  (C,)*“  -  k*  fl  i 

i= 1  i=  1 

where  C,  can  be  either  C"  or  C  bw;  Ra  is  defined  as  the  rate  of  the  aqueous  complexation  reaction. 
Adsorption/desorption  reactions  between  the  water  column  and  suspended  sediments:  n 6  [1,NJ,  ie 
[1,NC] 

C-W  +  s  **  C  s.  +  s 
s  _  usfc  p  w  _  u sb  q  r 

id  “  Kni  Kni  °iA 

where  Rads  is  defined  as  the  rate  of  the  adsorption  reaction. 

Adsorption/desorption  reactions  between  the  water  column  and  bed  sediments:  n6  [1,NJ,  i E  [1,NC] 

Ciw  +  Mn  *  C>Mn 

i1  - 

where  Radbl  is  defined  as  the  rate  of  the  adsorption  reaction. 

Adsorption/desorption  reactions  between  interstitial  water  and  bed  sediments:  n6  [1,NJ,  i6  [1,NC] 

Cjbw  +  Mn  -  Cnb  +  M„ 


2 
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where  Radb2  is  defined  as  the  rate  of  the  adsorption  reaction. 

Volatilization  reactions  between  the  water  column  and  atmosphere:  i£  [1,NC] 

Cw 

i  **  Pi 

R,  af.-,  w  ,  ab 

v  =  kj  Cj  -  kj  Pi 

where  Rv  is  defined  as  the  rate  of  the  volatilization  reaction.  In  Eqs.  (2.16),  (2.18),  (2.20),  (2.22), 
and  (2.24)  the  bold  form  is  used  to  present  the  chemical  formula  of  a  contaminant  or  a  sediment. 
To  achieve  transient  simulations,  the  concentrations  of  contaminants  and  sediments  must  be 
provided  initially.  Five  types  of  boundary  conditions  can  be  used  for  mobile  substances,  i.e., 
suspended  sediments,  dissolved  chemicals,  and  particulate  chemicals  on  suspended  sediments.  They 
are 

Dirichlet  boundary  condition: 

This  condition  is  applied  when  concentration  is  given  at  the  boundary.  That  is, 

C(xb,yb,t)  =  Cb(t) 

where  C  can  be  the  Sn,  C",  or  SnCnis;  Cb(t)  is  a  time-dependent  concentration  on  the  boundary  [  M/L’] . 
Variable  boundary  condition: 

This  boundary  condition  is  employed  when  the  flow  direction  would  change  with  time  during 
simulations.  Two  cases  are  considered,  regarding  to  the  flow  direction  on  the  boundary  segment. 
<  Case  1  >  Flow  is  coming  in  from  outside: 

b’t)  -  hK-V  C(xb,yb,t)J  = 

<  Case  2  >  Flow  is  going  out  from  inside: 

•  ( -  h  K  •  V  C(xb,yb,t) )  = 

where  n  is  the  outward  unit  normal  vector  of  the  boundary  segment;  Cin(t)  is  a  time-dependent 
concentration  [M/F3],  which  is  associated  with  the  incoming  flow. 

Ocean  boundary  condition: 
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This  boundary  condition  is  employed  when  the  flow  is  controlled  by  tides.  Two  cases  are 
considered,  regarding  to  the  flow  direction  on  the  boundary  segment. 

<  Case  1  >  Flow  is  coming  in  from  outside: 

/b,t)-hK-VC(xb,yb,t)j  = 

<  Case  2  >  Flow  is  going  out  from  inside: 

•  ( -  h  K  •  V  C(xb,yb,t) )  = 

where  CQ(t),  a  time-dependent  concentration  [M/L3]  at  the  ocean  boundary  during  the  incoming  flow, 
will  be  specified  as 

C  (t)  =  C  e  “ft 

where  Cmax  is  the  maximum  concentration  at  the  end  of  the  last  outgoing  flow  and  f  is  a  flushing 
decay  constant,  which  is  estimated  to  match  the  specific  flushing  rate  obtained  from  field  data. 
Cauchy  boundary  condition: 

This  condition  is  applied  when  the  total  incoming  material  flux  is  prescribed  as  a  function 
of  time  on  a  boundary  segment.  It  can  be  written  as 

where  qbc(t)  is  a  time-dependent  Cauchy  flux  [M/T/L]. 

Neumann  boundary  condition: 

This  condition  is  employed  when  the  incoming  diffusive  material  flux  can  be  prescribed  on 
a  boundary  segment.  It  is  written  as 

( -  hK- V  C(xb,yb,t))  =  < 


where  qbn(t)  is  a  time-dependent  Neumann  flux  [M/T/L]. 


Chapter  3.  NUMERICAL  APPROACHES 


In  this  chapter,  we  are  to  present  the  numerical  approaches  employed  to  solve  the  governing 
equations  of  flow  and  transport  given  in  the  previous  chapter. 


3.1.  Numerical  Approaches  for  Solving  the  Flow  Equations 

A  splitting  strategy  is  employed  in  our  approach  to  solve  2-D  depth-averaged  shallow  water 
equations.  After  mathematical  manipulation  Eqs.  (2.1)  through  (2.3)  can  be  written  in  matrix  form 
as 
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Eq.  (3.1)  can  be  written  as 
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With  the  Lagrangian  approach,  Eq.  (3.5)  can  be  further  written  as 


DE 

Dt 


RJ  +  RD 


where 


DE  dE  .  d  dE  .  d3E 
-  -  —  +  A  — +  A  — 

Dt  dt  dx  y  5y 


Eq.  (3.9)  can  be  approximated  as 


E  -  E  p 
At 


RJ  +  RD 


(3.6) 


(3.7) 


(3.8) 


(3.9) 


(3.10) 


(3.11) 


where  Ep  represents  the  Lagrangian  value  of  E.  If  we  let  E#  be  the  solution  of 
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we  have 


E#-Ep 

At 


then  substituting  Eq.  (3.13)  into  Eq.  (3.11)  yields 


E  -  E# 
At 


(3-13) 


(3.14) 


With  the  above  splitting  strategy,  we  solve  Eq.  (3.13)  for  E#  first  and  then  solve  Eq.  (3.14) 
to  obtain  E  based  on  E#.  We  also  employ  the  Picard  method  in  the  iteration  loop  that  deals  with  the 
coupling  of  Eqs.  (3.13)  and  (3.14)  due  to  non-linearity.  Eq.  (3.13)  is  solved  by  the  method  of 
characteristics  as  detailed  in  the  following. 

Method  of  Characteristics 

Eq.  (3.13)  can  be  written  as 


3E 

at 


+  A., 


3E  .  3E 

—  +  Av  — 

3x  y  3y 


RI 


(3.15) 


where 


RI, 

Ri 

RI  =  • 

ri2 

►  =  - 

R2  “  R^2 

Ri3 

r3  _  r^3 

(3.16) 


Let 
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Ak 


-  Axkx  +  Ayky 


ukx  +  vky  hkx  hky 

gkx  ukx  +  vky  0 

gky  0  ukx  +  vky 


(3.17) 


where  k  is  the  unit  vector  of  the  characteristic  direction.  Then  the  associated  eigenvalues  and 
eigenvectors  are 


*1  =  Ukx  +  Vky 

X2  =  ukx  +  vky+vgh 

X3  =  ukx  +  vky-Vgh 
We  define 


.  M. 

2  2  2 


■^h  gkx  gky 
2  2  2 


'  T 


The  inverse  of  matrix  L  is 


L  = 


L-1 


Tgh 

_  t/gh 

2 

2 

Skx 

gkx 

2 

2 

gky 

gky 

x  2 

2 

0 

ky 

-kx 

1 

K 

ky 

Vgh 

g 

g 

1 

K 

h. 

Vgh 

g 

g 

(3.18) 

(3.19) 


(3.20) 


(3.21) 


(3.22) 
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Define  the  following  transformation 


aw  =  L  _10E  = 


0 

1 


ky  -kx 
kx  ky 


t/gh 

1 


g  g 

K  K 

\Jgh  g  g 


0T1 

0U 

dv 


Multiplying  Eq.  (3.15)  by  L'1  yields 


aw  T  -i  A  T  aw  T  M  .  T  aw 

- +L  A  L - +  L  A.L - 


=  L  1 RI 


at 


dx 


3y 


or 


aw 

at 


u 

gcky 

2 

_  gcky 

2 

V 

_  §ckx 

2 

gckx 

2 

hky 

c 

u  +  ckx 

0 

aw 

+ 

dx 

_hkx 

c 

v  +  cky 

0 

_  hky 

0 

u  -ckx 

hkx 

0 

v-ck 

c 

c 

y 

aw 


ay 


=  L  _1  RI 


where 


c  =  Vgh 

Eq.  (3.25)  can  be  further  written  as  follows. 


(3.23) 


(3.24) 


(3.25) 


(3.26) 
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aw 

at 


+ 


u 

0 

0 

v  0 

0 

0 

u  +  ckx 

0 

dW 

+ 

0  v  +  ck 

0 

dx 

y 

0 

0 

u-ckx 

0  0 

v  -  cky 

aw 

ay 


\  ia_k  iV 

y  ax  x  dy 


+1  * 


k  k  ^_+k  k  _k  k 

y  y  ax  x  x  ay  x  y 


k  k  — +  k  k  —  -  k  k 

y  y  ax  x  x  ay  y 


} 

1  au  +  av  ^ 
V  ay  dxj 
( 


3u  +  dv 
dy  dx 


=  L  RI 


(3.27) 


Local  Diagonalization 

As  we  can  see  from  Eq.  (3.27)  which  displays  the  characteristic  form  of  the  shallow  water 
equations,  the  last  term  on  the  left  plays  the  role  of  a  non-linear  source/sink  and  may  increase 
difficulty  in  reaching  convergence.  The  method  we  are  taking  to  handle  this  is  to  make  the  term  zero 
using  diagonalization.  To  achieve  local  diagonalization,  we  first  define  a  new  L"1  which  performs 
the  following  transformation. 


dW  =  L _1 5E 


0 

c 


1 


dr] 

•  <3u  - 
dv 


egg 


(3.28) 


Multiplying  Eq.  (3.15)  by  this  new  L"1  will  yield 
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DWj 

ky(1)RI2-kx(1)RI3 

1  k(2)  ky(2) 

1ri1  +  ^ri2  +  ^ri3 

Dt 

DW, 

>  +  < 

s 

►  =  < 

Dt 

°2 

egg 

dw3 

,S3, 

- 1  k(2)  kf 

RIj  +  x  RI,  +  y  RI3 

.  c  g  g 

Dt 

where 


(3.29) 


s2  - 


k  (2)k  (2)  5u  +  k  (2)k  (2)  9V  _  k  (2) k  (2) 
y  y  dx  x  x  dy  x  y 


du  dv  ^ 

+ - 

dx  j 


3y 


S3  = 


i  (2),  (2)  da  ,  (2),  (2)  dv  (2),  (2) 

Ky  Ky  ~  +  Kx  Kx  —  Kx  Ky 

dx  dy 


du  dv 
dy  dx 


/  J 


(3.30) 

(3.31) 

(3.32) 

(3.33) 

(3.34) 

(3.35) 


It  is  clear  that  when  we  choose  (kx(1),  ky(1))  as  the  direction  for  the  first  characteristic  and  (kx<2), 
ky(2))  for  the  second  and  third  by  satisfying  the  following  conditions,  we  can  set  SI,  S2,  and  S3  to 
zero  and  thus  achieve  local  diagonalization. 


k 


(i)  dr)  =  0 

x  5y 


(3.36) 
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,  (2),  (2)  du  ,  (2),  (2)  dv  _  v  (2),  (2) 

Ky  Ky  -  Kx  Kx  _  Kx  Ky 

dx  dy 


/  du  +  dv  N 
dx 


dy 


=  0 


(3.37) 


We  can  always  find  kx(1)  and  ky(1)  to  satisfy  Eq.  (3.36),  i.e.,  wave  direction  that  is  tangential 
to  the  gradient  of  water  surface  elevation  can  always  be  found,  but  we  will  not  be  able  to  find  any 
real  solutions  of  kx(2)  and  ky(2)  in  solving  Eq.  (3.37)  when  the  discriminant  of  the  quadratic  equation, 
i.e.,  (du/dy  +  <3v/<3x)2  -  4(<3u/<3x)(<3v/<3y),  is  less  than  zero.  In  this  case,  Hirsh  et  al.  have  suggested 
to  select  directions  which  minimize  the  coupling  terms  [Hirsh  et  al,  1987],  e.g.,  S2  and  S3  in  our  case 
here.  Therefore,  we  calculate  kx(2)  andky(2)  through  the  following  equation  [Hirsh  et  al,  1987]  if  there 
exists  no  real  (kx(2),  ky(2))  satisfying  Eq.  (3.37). 


tan 


2k 


(2) 


.  (2) 


du  dv 
dy  dx 
du  _  dv 
dx  dy 


(3.38) 


Backward  Tracking  Along  Characteristics 

Eqs.  (3.30)  through  (3.32)  can  be  discretized,  in  time,  as  follows. 


,(i)Un+1-<  k(Dvn+1-v 


At 


At 


-  +  Sj  =  k(1)RI2-kx(1)RI3 


(3.39) 


n  + 1  *  i  (2)  n + 1  *  i  (2)  n  + 1  * 

1  n  -n2  kx  U  “U2  ky  V  “V2  - 

- + - +  — i - +  So 


.(2) 


.(2) 


c  At  g  At  g  At 
_  \  nn+1  -ti*  kx(2)  un+1  -u3*  ky(2)  vn+1  -v3* 


At 


g  At  g  At 


+  So  = 


-  RIj  +  —  RI2  +  RI3  (3 .40) 

eg  g 


(2)  (2) 

—  RIj  +— RI2  +  -^-RI3  (3.41) 

c  g  g 


where  u,  and  v , '  are  determined  through  backward  tracking  along  the  first  characteristic,  while  r)2*, 
u2*,  and  v2*  are  associated  with  the  second  characteristic  and  T)  (  ,  u3*,  and  v3*  are  associated  with  the 
third  characteristic.  Figure  3.1  demonstrates  this  backward  tracking  along  characteristics.  In 
general,  we  have 
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* 

U1 

=  aiUkl 

+  a2Uk2 

+  a3  Uk3 

+  a4Uk4 

(3.42) 

* 

V1 

=  aiVkl 

+  a2Vk2 

+  a3Vk3 

+  a4Vk4 

(3.43) 

* 

*1  2 

1  n 
=  b.-Hj, 

+b2np 

+  b3  *1j3 

+  b4  *lj4 

(3.44) 

* 

U2 

=  blUjl 

+  b2Uj2 

+  b3Uj3 

+  b4Uj4 

(3.45) 

* 

V2 

=  bivii 

+  Vj2 

+  Vj3 

+  V; 

(3.46) 

* 

*13  = 

dl*lml  +  d2*1m2 

+  d3  *1  m3  +  d4  *1  m4 

(3.47) 

* 

u3  = 

i  n 

;  dlUml 

+  d2Um2 

J  11 

+  d3Um: 

1  +  d4Um4 

(3.48) 

* 

V3  = 

i  n 

:  Vml 

+  d2Vm2 

J  11 

+  d3Vm3 

!+d4Vm4 

(3.49) 

where  Tl/1,  u",  v,"  represent  the  water  surface  elevation,  x-velocity,  and  y-velocity  at  node  i  (i  =  kl 
through  k4,  jl  through  j4,  and  ml  through  m4)  at  t  =  tn  (Figure  3.1);  coefficients  a,  through  a4  are 
interpolation  factors  at  nodes  kl  through  k4,  respectively,  while  h,  through  b4  are  corresponding  to 
nodes  j  1  through  j4  and  d,  through  d4  are  for  nodes  ml  through  m4. 
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jl  j2 

Figure  3.1.  Backward  tracking  along  characteristics  in  solving 
2-D  hydrodynamic  equations. 

By  using  the  spatial  reachout  scheme  [Lai,  1987]  in  backward  tracking,  u,  \  v*,  r)2*,  u2*,  v2*, 
r|3*,  u3*,  and  v3*  can  be  determined  explicitly  as  described  by  Eqs.  (3.42)  through  (3.49).  Hence,  we 
can  solve  Eqs.  (3.39)  through  (3.41)  for  (T)#,  u#,  v#)  node  by  node  as  stated  below. 

After  being  multiplied  by  water  depth  (i.e.,  h),  Eqs.  (3.39)  through  (3.41)  can  be  written  as 


ai2U  +  ai3V  =bi 

(3.50) 

+a22U+a23V=b2 

(3.51) 

a3l  p  +a32u  +  a33v  =  b3 

(3.52) 

where 
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ai2 


(1) 


—  +hK  +  S 
At 


(3.53) 


a!3  =  -kx 


(1) 


— +hK+S 
At 


(3.54) 


bi  = 


hk 


(i) 


hk 


At 


-ui 


(i) 

x  * 


At 


V!  -gh 


/  ^  (i)  _  k  (i)  3r)  N 


dx 


dy 


+  k(l)u  sS  -  k(1|v  sS  +  k<l!—  -k<1!— 
y  x  y  _  x 


(3.55) 


21 


c  At 


(3.56) 


a22  kx 


(2) 


h  +i(hK+S) 


a23  k 


(2) 


gAt  g 


h  +i(hK  +  S) 


bo  = 


hk 


(2) 


hk 


gAt  g 

(2) 


h  dZQ 
c  dx 

h  dZp 
c  dy 


— n2  +  — 

cAt  gAt 


U2  + 


,  k<2)  v(2) 

— —  v2*  +  ~  SS  +  — u  s  S  +  — v  s  S 
gAt  c  g  g 


k(2)  s  (2)  s 

Tx  [  Ky  Ty  h 

g  P  g  P  c 


,  (2),  (2)  du  ,  (2),  (2)  dv  _  ,  (2),  (2)| 
y  y  dx  x  x  dy  x  y 


5u  +  dv 
dy  dx 


)  J 


(3.57) 


(3.58) 


(3.59) 


a31 


-h 
c  At 


a32  =  k 


(2) 


h  +  — (hK  +  S) 


gAt  g 


h  5zo 
c  dx 


(3.60) 

(3.61) 


a33  =  k 


(2) 


h  +  — (hK  +  S) 


gAt  g 


h  dZ0 
c  3y 


(3.62) 
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b,  = 


-h 
c  At 


hk 


(2) 


hk 


(2) 


gAt 


u3  + 


h  k 


(2) 


k(2)  ‘  k(2)  » 

Kx  Tx  Ky  Ty 
+ - +  - 3-  + 

g  P  g  P 


gAt 
h2 


v3  -  ~S  + 


u  + 


c(2) 

—  v  sS 


(3.63) 


,  (2),  (2)  du  ,  (2),  (2)  5v 
y  y  dx  x  x  dy 


k(2)k(2) 

Kx  Ky 


/  5u  +  3v  ^ 
v  3y  3xy 


With  the  coefficients,  arj,  and  the  right  hand  side,  b,.  evaluated  using  the  updated  iterate,  we 
solve  Eqs.  (3.50)  through  (3.52)  for  T)#,  u#,  v#  arithmetically. 

Finite  Element  Approximation 

Now,  we  can  solve  Eq.  (3.14)  as  follows.  Eq.  (3.14)  can  be  written  in  detail  as 


At 


(3.64) 


h(u-u#) 

At 


dF  dF 

xx  +  yx 

dx  dy 


(3.65) 


h(v  -  v#) 
At 


3Fxy  dF 


yy 


dx  dy 


(3.66) 


Eq.  (3.64)  yields  that  T)  =  T|#,  while  Eqs.  (3.65)  and  (3.66),  in  the  form  of  diffusive  transport 
equations,  can  be  discretized  with  the  Galerkin  finite  element  method  as 

[QA]  {u}  +  [QD  1  ] { u }  +  [QD4]  { v}  =  {QR’J  +  IBQ1}  (3.67) 


— —  [QA]  { v}  +  [QD2]  {u}  +  [QD3]  {v}  =  {QR2}  +  {BQ2} 
At 


(3.68) 


where  [QA]  is  the  mass  matrix;  [QD1]  and  [QD3]  are  stiffness  matrices  involving  only  the  diagonal 
terms  of  eddy  diffusion;  [QD2]  and  [QD4]  are  stiffness  matrices  involving  only  the  off-diagonal 
terms  of  eddy  diffusion;  [QR1 }  and  [QR2]  are  load  vectors  generated  from  u#  and  v#;  [BQ1 }  and 
[BQ2]  are  load  vectors  representing  diffusive  fluxes  through  the  boundary.  They  can  be  written  in 
detail  as 
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QAij  =  E  f  N„e  h  Np6  dR„ 

eeM„  " 


v-  ,  c  he**  o 

QDij  =  E  f  vn;- 

eeM„  ^ 


0  he 


xy 


VNp  dRe 


QD,j  =  E  fVN.-- 

eeM„  " 


0  he 
0  0 


xy 


VNp  dRe 


QD,j  =  E  f  VN„e- 

eeM„  i 


h€xy  0 

0  he 


yy 


V  Np  dRe 


eeM„ 


/vn;- 

0 

o' 

he 

0 

R„ 

L  xy 

VNp  dRe 


— —  dRe 

eeM  £  At 


QR,1  =  E  f  N.' 

PCM  J 


QRi2  =  E  f 

ocA/f 


N."  ^  dRs 


eeM„  ^  At 

R„ 


( 

he  0 

BQ,1  =  E  f  N„V 

e6Me  t  1 

De 

XX 

.  0  h6xy 

V  u  +  N0  n 


0  0 
h6xy  0 


Vv 


r  dB„ 


BQ2  =  E  {■ 

eeMp  i 


xt  e 

Nan 


0  he 
0  0 


xy 


V  u  +  N0  n 


h6xy  0 

0  he 


yy 


Vv 


[  dB„ 


(3.69) 

(3.70) 

(3.71) 

(3.72) 

(3.73) 

(3.74) 

(3.75) 

(3.76) 

(3.77) 


where  Me  represents  all  elements  composing  the  domain  of  interest;  Re  represent  the  domain  of 
element  e;  Be  represents  the  boundary  of  element  e. 

Implementation  of  Boundary  Conditions 
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A  boundary  segment  can  be  either  a  closed  (impermeable)  or  an  open  (flow-through) 
boundary.  To  perform  numerical  simulations,  no  boundary  condition  needs  to  be  specified  on  the 
closed  boundary,  while  T),  u,  and  v  of  the  approaching  flow,  with  respect  to  the  characteristic 
ongoing  direction,  must  be  given  on  the  open  boundary.  In  solving  Eq.  (3.13),  we  implement 
boundary  conditions  after  using  backward  tracking  to  construct  Eqs.  (3.39)  through  (3.41)  as  follows. 
For  the  case  of  a  closed  boundary  node: 

The  first  characteristic  equation  is  replaced  by  the  following  equation,  representing  a  zero 
normal  flux  on  the  boundary. 

n-V  =  0  (3.78) 


Tracking  along  the  impermeable  boundary  will  be  used  for  the  second  and  third  characteristics  if 
necessary. 

For  the  case  of  a  open  boundary  node: 

In  this  case  we  assume  a  zero  lateral  flux  on  the  boundary,  and  the  first  characteristic 
equation  is  replaced  by  the  following  equation. 

n±V  =  0  (3.79) 

where  n  is  the  unit  vector  parallel  to  the  boundary.  Either  the  second  or  the  third  characteristic 
equation,  depending  on  the  propagation  direction,  will  be  replaced  by  the  boundary  condition 
specified  on  the  boundary  where  sub-critical  flow  is  presumed.  Only  when  the  propagation  direction 
is  inward  to  the  domain  of  interest,  will  the  replacement  be  employed.  There  are  two  types  of 
boundary  conditions  considered  here. 

(1)  Given  time-dependent  normal  flux,  qb(t,  xb,  yb): 

n-V  =  —  (3.80) 

h 

(2)  Given  time-dependent  water  surface  elevation,  T)b(t,  xb,  yb): 


n  = 


(3.81) 
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Based  on  the  fact  that  diffusive  fluxes  are  negligible  when  compared  with  advective  fluxes  in 
bay/estuary  areas,  we  have  assumed  zero  diffusive  fluxes  through  open  boundaries  in  using  the 
Galerkin  finite  element  method  to  solve  Eq.  (3.14).  This  assumption  will  lead  to  both  {BQ1 }  and 
{ BQ2 }  being  zero,  which  consequently  eases  the  construction  of  the  matrix  equations  corresponding 
to  Eqs.  (3.65)  and  (3.66). 


3.2.  Numerical  Approaches  for  Solving  the  Transport  Equations 

In  bay/estuary  areas  where  advection  dominates  dispersion,  the  LE  approach  [Wood  and 
Baptista,  1992]  has  been  considered  an  appropriate  way  to  compute  material  transport  because  it  may 
greatly  reduce  many  types  of  numerical  errors  [Yeh,  1990] .  This  approach  can  always  provide  non¬ 
negative  computed  concentration  in  solving  solute  transport  equations  when  linear  or  bi-linear 
elements  are  used  for  discretization,  which  is  required  in  computing  concentration  distributions 
subject  to  chemical  reactions  [Cheng  and  Yeh,  1998].  Therefore  we  can  write  the  governing 
equations  of  mobile  substances  in  the  LE  form  as  follows. 

For  suspended  sediments:  n6  [1,NJ 


dS 


dt 


n  -  VjhK-VS„ 


Mn+Rn-Dn 


-  ss„ 


(3.82) 


where  d/dt  represents  the  material  or  total  derivative  with  respect  to  time;  S  is  external  fluid  source 
[L/T], 

For  dissolved  chemicals  in  the  water  column:  i6  [1,NC] 
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dCj 

h — — -V- 
dt 


hK-VC; 


rx 

E  (ami  -bmi)kmh 


m=l 


Nc  ,  rf  Nc 

n  <cDb--^n  (Cj*)' 


mj 


(3.83) 


R 


M  ™  +  hk  *bp,  +  £  k „?h (S„  C„’)  +  £  kbh  (MnC„b)  +  M,"  +  Eq1”  +  £  — 0,0, 

n  Pn 


bw 


n= 1 


n  =  l 


D. 


A'h+hki,,+x;s]1hkj+x;Mnk1b'+E+5]^0n+s 

n  Pn 


n  =  l 


n  =  l 


c, 


For  particulate  chemicals  on  suspended  sediments:  nE  [1,NJ,  iE  [1,NC] 


dt 


hK-V(SnO 


K + knfhsnqw  +  ^(MnO  -  -^csnO 


R_ 


D„ 


M 


(3.84) 


0+k>  +  s 


(SnO 


In  obtaining  Eqs.  (3.82)  through  (3.84),  we  have  substituted  the  flow  continuity  equation  into 
Eqs.  (2.10),  (2.12),  and  (2.13).  The  flow  continuity  equation  is  the  balance  equation  of  water  mass 
and  can  be  written  as  follows. 


3ri  3uh  3vh 

- L  + - + - 

3t  dx  dy 


(3.85) 


where  T)  =  T|(x,y,t)  is  water  surface  elevation;  u  is  the  x-velocity  [L/T];  v  is  the  y-velocity  [L/T]; 
Here  we  also  define  sea  bed  elevation  Z0  =  -d(x,y),  where  d(x,y)  is  water  depth  at  mean  sea  level. 
It  is  thus  obvious  that  ri(x,y,t)  =  Z0(x,y)  +  h(x,y,t)  and  dr)/<3t  =  dh/<3t. 

In  solving  Eqs.  (3.82)  through  (3.84),  we  compute  advection  in  the  Lagrangian  step  where 
backward  particle  tracking  is  applied  to  determine  the  Lagrangian  concentrations  at  all  global  nodes. 
We  then  use  the  Lagrangian  concentrations  as  the  initial  conditions  to  solve  the  transport  equations 
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in  the  Eulerian  step  where  the  Galerkin  finite  element  is  applied. 

The  2-D  transport  equations  presented  above  are  a  set  of  non-linear  partial  differential 
equations  where  non-linearity  is  introduced  through  the  deposition/erosion  and  reaction  terms.  To 
perform  efficient  computation  without  sacrificing  much  computational  accuracy,  we  employ  a 
predictor-corrector  strategy  to  solve  the  governing  equations.  Here  we  give  a  brief  description  of 
this  strategy.  The  governing  equations  of  mobile  substances  (i.e.,  Eqs.  (3.82)  through  (3.84))  can 
be  expressed  in  the  following  form. 

h~  =  L(C)  +  RHS  (3.86) 

dt 


where  C  can  be  Sn,  Qw,  or  SnCnis,  L  is  the  dispersion  operator,  and  RHS  represents  other  terms.  With 
the  predictor-corrector  strategy,  we  solve  Eq.  (3.86)  with  the  following  two  steps.  First,  we  solve 

riNtl/2_p  * 

h— - —  =  L(Cn  +  1/2)  +  (RHS)n  (3.87) 

Ax 


where  Ax  represents  the  time  period  consumed  in  backward  particle  tracking  [T];  C  is  the 
Lagrangian  concentration  [M/L3] ;  (RHS)N  represents  RHS  evaluated  at  the  previous  time  [M/L2/T] ; 
CN+1/2  is  the  intermediate  concentration  which  we  solve  for  with  RHS  evaluated  at  the  previous  time 
[M/L3] .  It  is  noted  that  At  is  the  real  time  period  for  a  fictitious  particle  to  travel  through  advection 
from  a  upstream  location  to  the  global  node  being  considered,  Ax  is  location-dependent  [Yeh  et  al, 
1992].  This  time  period  is  also  what  we  use  to  calculate  dispersion  and  reactions  in  the  Eulerian 
step.  Second,  we  solve 

hCN+1~C  =  L(Cn  +  1/2)  +  (RHS)n+1  (3.88) 

Ax 


where  CN+1  is  the  concentration  of  the  present  time  [M/L3];  (RHS)N+1  represents  RHS  evaluated  at 
the  present  time  [M/L2/T].  We  can  subtract  Eq.  (3.88)  by  Eq.  (3.87)  to  yield 


pN  +  l 

h— - 


-  N  +  1/2 


Ax 


=  (RHS)n  +  1  -(RHS)n 


(3.89) 
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To  save  computer  effort,  it  is  better  to  solve  Eq.  (3.89)  than  Eq.  (3.88)  because  Eq.  (3.89) 
does  not  involve  the  dispersion  operator  and  thus  can  be  solved  node  by  node.  To  summarize,  we 
first  solve  Eq.  (30)  for  CN+1/2  in  the  predictor  step  and  perform  the  corrector  step  by  solving  Eq. 
(3.89)  to  obtain  CN+1. 

We  employ  an  adaptive  explicit-implicit  scheme  to  avoid  negative  numerical  results  in 
solving  Eq.  (3.87).  With  this  scheme,  the  following  equation  is  actually  solved  when  (RHS)N  is 
negative. 


p N  + 1/2  n  *  rr?rjQs\N 

frW  ^N+l/2  _  L('Cn  +  1/2') 

'  At  c  n 


(3.90) 


Consequently,  we  are  solving  the  following  correction  equation  in  the  corrector  step. 

qN+1  _  ^N  +  l/2 


At 


(RHS) 


N  + 1 


(RHS)n 


iN  +  1/2 


(3.91) 


With  the  Galerkin  finite  element  approximation,  Eq.  (3.87)  can  be  discretized  as  follows. 


-^i  +  w  [QD] 
At  j 


{C}  =  I^1{CM  +  {QR}+(1-w)[QD]{Cn}  +  {QB}  (3.92) 
At 


where  w  is  time  weighting  factor;  [QA]  is  the  mass  matrix,  [QD]  is  the  stiffness  matrix  due  to 
dispersion;  {QR}  is  the  load  vector  contributed  from  (RHS)N;  { QB }  is  the  load  vector  accounting 
for  dispersive  boundary  fluxes.  The  associated  element  matrices  and  column  vectors  are  defined  as 
follows. 


[QA]e  =  |NihNjdRe 

Re 

[QD]e  =  JvNi-(hK-VNj)dRe 

Re 

{QR}e  =  j‘Ni(RHS)NdRe 

Re 


(3.93) 

(3.94) 

(3.95) 
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{ QB  }e  -  -  J(-hK-VC)dBe 
Be 


(3.96) 


where  Re  and  Be  represent  the  domain  and  the  boundary,  respectively,  of  element  e;  N,  and  Nj  are 
element  shape  functions  at  the  i-th  and  j-th  element  nodes,  respectively,  of  element  e.  It  is  noted  that 
{QR}e  must  be  evaluated  using  nodal  quadratures  such  that  Eqs.  (3.87)  and  (3.89)  are  solved 
consistently.  Also,  boundary  conditions  for  mobile  substances  are  implemented  in  this  predictor 
step.  If  (RHS)n  is  negative,  Eq.  (3.90)  is  solved  in  the  predictor  step.  It  can  be  discretized  as 

^!  +  w[QD]-[QC]  )  {C}  =  ^  { C  * }  +  ( 1  -  w )  [  QD  ]  { C  N }  +  {  QB  }  (3.97) 

(Ax  j  Ax 


where  [QC]  is  the  stiffness  matrix  due  to  source/sink  terms  evaluated  at  the  previous  time.  The 
element  matrix  of  [QC]  is  defined  as 


[QC]'  =  |n, 


<™^NjdRe 

ri  N  1 


(3.98) 


where  [QC]e  must  be  lumped  to  maintain  consistency  in  solving  Eqs.  (3.90)  and  (3.91). 

With  the  predictor-corrector  strategy,  we  solve  the  entire  transport  system  by  completing  the 
following  steps. 

Step  1  We  determine  the  Lagrangian  concentrations  for  suspended  sediments,  dissolved  chemicals 
in  the  water  column,  and  particulate  chemicals  on  suspend  sediments,  respectively. 

Step  2  We  solve  Eqs.  (3.82)  and  (2.1 1)  for  Sn  and  Mn,  respectively.  In  this  step,  we  first  solve  Eq. 
(3.82)  with  all  source/sink  terms  evaluated  at  the  previous  time  to  obtain  the  intermediate  value  for 
suspended  sediments.  Then  we  prepare  the  corrector  form  of  Eq.  (3.82)  as  follows. 
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h  (Sn)N+1~(Sn)N+1/2 

Ax 

=  (RHSns)N+1 -(RHSns)N 

s  M  ,  (RHSns)N 
-  (RHSn  )N+1  - - ~  ( Sn ) 

(Sn) 

where 

(RHSns)  =  [Mns  +  Rn-Dn]-SSn  (3.100) 

At  last,  we  solve  Eqs.  (3.99)  and  (2)  node  by  node  with  the  Picard  method  to  obtain  the  solutions  for 
suspended  and  bed  sediments. 

Step  3  We  solve  Eqs.  (3.83),  (3.84),  (2.14),  and  (2.15)  forQw,  SnCnis,  MnCm\  andCnibw,  respectively. 
In  this  step,  we  solve  Eqs.  (3.83)  and  (3.84)  with  all  source/sink  terms  evaluated  at  the  previous  time 
to  determine  the  intermediate  concentrations  of  C"  and  SnCnis,  followed  by  preparing  the  corrector 
form  of  Eqs.  (3.83)  and  (3.84)  as  given  below. 

For  dissolved  chemicals: 


(C-W)N+1  -  (CW)N+1/2 


Ax 

=  (RHS;C)n+1  -(RHSiC)N 

(if  (RHSiC)N>o) 

(3.101) 

=  (RHSjC)n+1  - (RHSi  )N  (C;W)N+1/2 

(if  (RHS.°)n<0  ) 

\  1  / 

(C;W)N 

where 
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N_ 


(RHS,C)  =  E  <am|-bmpk*h 


m=l 


Nc  k  rf  N( 

n  (c/)b--^n  (Cjw)’ 

j=i 


i  rb  .  1 
km  J  =  1 


N. 


N. 


N. 


R 


M  ”  +  hk,,lPi  +  E  k>(S„  C„*)  +  E  C  (MnC  j)  +Mi”  +  ECi6w  +  E  — 1 0„Ci 

n  =  l  Pn 


bw 


n  =  l 


n  =  l 


(3.102) 


N. 


D„ 


X'h  thk-'+E  Snhk^'  +E  Mnknif+S+E  +  E  — 8„ 

n  =  l  Pn 


n  =  l 


n  =  l 


c; 


For  particulate  chemicals  on  suspended  sediments: 


/  q  s  \N  + 1  _  /  q  /~*  s  \N  +  1/2 

h - 


At 

=  (RHSnj)N+1  -  (RHSnS;)N  ( 

if  (RHSn;)N^o) 

(3.103) 

_  /puc^N+l  (RHSni)  s.n+1/2 

(if(RHSn*)N<o) 

(O 


where 


(RHSn*)  - 


M“  +  knfhSnC,"  +  -^-(MnC„1;)  -  ^(S„C„5 
M„  S„ 


^h  +  knSiDh  +  S 


(SnO 


(3.104) 


Eqs.  (2.14)  and  (2.15)  can  be  written  as  follows. 
For  particulate  chemicals  on  bed  sediments: 


(MnC^)N+1-(MnCn^)N 

At 


(RHSV*1 


(3.105) 


where 
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(RHSn  )  - 


+  k  "m„C,'  +  k>„C, 


bw 


M 


,  b  ,  bb  ,  2b 

A ■ +k  ■  +k  • 

m  m  m 


(MnCn") 


(3.106) 


For  dissolved  chemicals  in  interstitial  water  of  bed  sediments: 


Nc 


0. 


\  N  +  l  /  N  \  N 


E-M.C, 

V  n=1  Pn 


bw 


E— MnC, 

\ n=1  Pn 


bw 


At 


(RHS;bw)N+1 


(3.107) 


where 


N_ 


(RHSjbw)  =  E  k>nC„-  - 


n= 1 


Ns 


D„ 


E  +  E— 0,  +  Ekni'MB 

y  n=l  Pn  n=l 


C,bW  + 


Ns  D  X 

e  +  E— 0.Q' 

v  n=l  Pn 


N, 


0 


E^m„ 

n  =  l  Pn 


N„ 


m=l 


E  (amj-bmj)C 


Nc  ,rf  Nc 

n  (Cjbw)b-^n  (c;w)a 

j^1  Ch 


(3.108) 


Finally,  we  solve  Eqs.  (3.101),  (3.103),  (3.105),  and  (3. 107)  node  by  node  with  the  Newton-Raphson 
method.  When  the  chemistry  system  is  highly  non-linear,  we  may  have  Jacobian  entries  vary  many 
order  of  magnitude,  which  will  introduce  numerical  difficulty  in  reaching  correct  solutions  [Press 
et  al.,  1992].  To  overcome  this,  we  incorporate  full  pivoting  with  a  direct  solver  to  correctly  solve 
the  Jacobian  matrix  equation  resulting  from  the  Newton-Raphson  method  [Yeh  et  al.,  1998].  The 
details  of  the  residual  equations  and  the  associated  Jacobian  matrix  for  the  Newton-Raphson  method 
are  given  below. 

To  achieve  Step  3  above,  the  following  residual  equations  are  constructed,  where  we  have 
i  6  [1,  Nc]  and  n  E  [1,  NJ  for  the  ranges  of  subscripts  i  and  n. 

For  dissolved  chemicals  in  the  water  column: 
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(RES,0) 


=  h 


(C-W)N+1  -  (Cw)N  +  ,/2 


-  At 


(RHSjC)n+1  -  (RHSjT  (if  (RHSjT  *  0 


<=\N 


CxN. 


=  h 


W\N  + 1  //-i  W\N  +  1/2 


(q  r+l  -(Cj  > 


■At 


C\N 


(RHSj  ) 


+  1  (RHSjT 


WxN 


(Ci  > 


w  \N  + 1/2 


(Cj")- 


if  (RHSiC)N<0 


For  particulate  chemicals  on  suspended  sediments: 


(RES®) 


=  h 


(sncn-)N+1  -(sncn;> 


S\N  +  1/2 


■  At 


if  (RHSni)N>0 


(RHSni)N+1-(RHSni)N 


=  h 


(SnCiNtl-(SnQ 


S\N  +  1/2 


■  At 


(RHSn) 


SxN+ 


,  (RHSn*)N 


s  \N 


(Sncn;) 


SxN  +  1/2 


if  (RHSni)N<0 


(SnCni> 


For  particulate  chemicals  on  bed  sediments: 


(RES^)  =  (MnC^)N+1-(MnC^)N-At(RHS^)N+1 


For  dissolved  chemicals  in  the  interstitial  water  of  the  bed  sediments: 


(RESjbw)  = 


'  Ns  AN+1 

E—  M„Cib' 
v  n=1  Pn  ) 


N. 


0, 


E  — Mn  cr 

n=l  Pn  j 


N 


■  At  (RHSibw)N+1 


The  entries  of  the  associated  Jacobian  matrix  are  evaluated  as  follows. 


(3.109) 


(3.110) 


(3.111) 


(3.112) 


For  dissolved  chemicals  in  the  water  column: 
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(RES;0) 


:cr> 
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(3.115) 

(3.116) 


For  particulate  chemicals  on  suspended  sediments: 
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a(REsn;) 
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(3.117) 


(3.118) 


(3.119) 


(3.120) 


For  particulate  chemicals  on  bed  sediments: 
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(3.122) 


(3.123) 


(3.124) 


For  dissolved  chemicals  in  the  interstitial  water  of  the  bed  sediments: 
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(3.128) 


3.3  Estimation  of  Deposition  and  Erosion 

We  estimate  the  deposition  and  erosion  rates  by  using  the  following  equations. 
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<  Case  1  >  Cohesive  sediments:  (e.g.,  silt  and  clay) 


D 


n 


=  Y  S 

sn  n 


(3.129) 


R 


n 


(3.130) 


where  Vsn  is  the  settling  velocity  of  the  n-th  size  fraction  sediment  [L/T] ;  Tb  is  the  bottom  shear  stress 
or  the  bottom  friction  stress  [M/L/T2];  TcDn  is  the  critical  shear  stress  for  the  deposition  of  the  n-th 
size  fraction  sediment  [M/L/T2];  XcRn  is  the  critical  shear  stress  for  the  erosion  of  the  n-th  size 
fraction  sediment  [M/L/T2] ;  En  is  the  erodibility  of  the  n-th  size  fraction  sediment  [M/L2] .  In  the 
computer  code,  Vsn,  XcDn,  XcRn,  and  En  are  input  parameters,  while  Th  is  computed  in  the  flow  module. 
<  Case  2  >  Non-cohesive  sediments:  (e.g.,  sand) 


G4  -G 

q  _  sAn  sn 

n  AL 


G  -G  A 

_  sn  sAn 

n  AL 


(3.131) 


(3.132) 


where  GsAn  is  the  actual  load  rate  of  the  n-th  size  fraction  sediment  per  unit  width  at  a  upstream 
location  [M/L/T] ;  Gsn  is  the  maximum  load  rate  (capacity)  of  the  n-th  size  fraction  sediment  per  unit 
width  at  a  downstream  location  [M/L/T];  AL  is  the  distance  between  the  upstream  and  the 
downstream  locations.  In  the  computer  code,  AL  can  be  determined  based  on  the  coordinates,  while 
GsAn  and  Gsn  are  computed  based  on  the  following  equations. 

GsAn  =  Sn*u*r  (3.133) 


p2urS(xb-xcm) 

gdn(Psn-p)2 


(3.134) 


where  p  is  the  fluid  density  [M/L3];  psn  is  the  density  of  the  n-th  size  fraction  sediment  [M/L3];  u  is 
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channel  flow  velocity  [L/T] ;  r  is  hydraulic  radius  [L] ;  S  the  friction  slope,  dn  is  the  median  diameter 
of  the  n-th  size  fraction  sediment  particle  [L] ;  g  is  gravity  [L/T2];  Tcrn  is  the  critical  bottom  shear 
stress  of  the  n-th  size  fraction  sediment  at  which  sediment  movement  begins  [M/L/T2] .  Among  these 
parameters,  p,  psn,  dn,  and  Tcrn  are  input  by  users,  while  u,  r,  and  S  are  estimated  in  the  flow  module. 
Note:  Eq.  (3.134)  is  from  “Hydraulics  of  Sediment  Transport”  by  Walter  Hans  Graf  (1984),  Water 
Resources  Publication,  Eq.  (7.14)  on  page  128.  Eqs.  (3.129)  through  (3.132)  can  be  found 
from  “CHNTRN:  A  Channel  Transport  Model  for  Simulating  Sediment  and  Chemical 
Distribution  in  a  Stream/River  Network”  by  Yeh  (1983),  ORNL-5882  Report,  Eqs.  (29)  and 
(30)  on  page  12. 
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Chapter  4.  EXAMPLES 

4.1  Example  1  —  1-D  Standing  Wave  Problem 

We  used  a  one-dimensional  standing  wave  example  [Wang  and  Connor,  1975]  to  verify  our 
approach.  This  problem  was  governed  by  the  wave  equations  as  follows. 

02u  _  c2  S2u 

0  t 2  0  X  2 

d2i)  =  c2  52T1 
0  t 2  0  X  2 

c2  =  gh 

where  u  is  the  velocity  along  the  flow  direction  [L/T];  c  is  the  wave  speed  [L/T];  h  is  water  depth 
[L] ;  T|  is  the  deviation  of  water  surface  elevation  [L] ;  g  is  gravity  [L/T2] .  The  domain  of  interest  was 
200  m  long  in  the  x-direction  and  50  m  wide  in  the  y-direction.  It  was  discretized  with  20  elements: 
10  m  x  50  m  each.  The  mean  water  depth  was  4  meters.  The  boundary  end  at  X  =  200  m  was 
closed,  while  the  other  one  at  X  =  0  m  was  open  and  the  water  surface  elevation  fluctuated  up  and 
down  according  to 

,  .  2  7C  t 


where  T|0  represents  the  amplitude  of  the  sinusoidal  surface  deviation  applied  to  the  boundary  X  = 
0  m.  Therefore,  we  set  u  =  0  m/s  at  X  =  200  m  to  settle  the  boundary  condition  for  achieving  the 


simulation.  The  analytical  solution  of  this  linear  wave  problem  is  as  follows. 
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Given  initial  condition  at  t  =  50  s  according  to  Eqs.  (4.5)  and  (4.6),  we  first  chose  T|0  =  0.1 
in  this  example  to  perform  a  simulation  of  400  seconds,  i.e.,  from  t  =  50  s  to  450  s.  We  found  the 
numerical  solution  is  essentially  unchanged  for  time  step  size  less  than  2.5  seconds  and  element 
number  more  than  20,  which  means  we  have  obtained  a  convergent  solution  to  the  hydrodynamic 
flow  equations  (i.e.,  Eqs.  (2.1)  through  (2.3)).  The  differences  between  the  numerical  result  and  the 
analytical  solution  of  Eqs.  (4.1)  through  (4.3)  (  (Figures  4.1(a)  and  4.2(a))  was  caused  by  the 
nonlinear  terms  accounted  for  in  Eq.  (2.2).  To  prove  this,  we  reduced  T|0  to  0.05  and  0.025  to 
diminish  the  nonlinear  effect.  Figures  4. 1  and  4.2  compare  the  numerical  results  and  the  analytical 
solutions  of  the  deviation  of  water  surface  elevation  at  x  =  100  m  and  200  m,  respectively  for  r|0  = 
0.1,  0.05,  and  0.025.  It  is  obvious  that  the  difference  decreases  with  the  reduction  of  T|0.  When  T|0 
is  reduced  to  0.025,  we  have  the  computational  result  and  the  analytical  solution  match  precisely 
(Figures  3(c)  and  4(c)).  An  excellent  agreement  was  also  obtained  in  comparison  of  the  velocity  at 
time  =  100  s,  200  s,  300  s,  and  400  s  for  the  case  of  T|0  =  0.025  (Figure  4.3).  It  is  thus  verified  that 
our  MOC  numerical  model  can  solve  standing  wave  problems  accurately. 


Figure  4.1.  The  deviation  of  water  surface  elevation  at  x  =  100  for  (a)  T|0  =  0. 
(b)  T|0  =  0.05  m,  and  (c)  r|0  =  0.025  m  during  the  simulation  of  Example  1 


(b)  T|0  =  0.05  m,  and  (c)  T|0  =  0.025  m  during  the  simulation  of  Example  1 


x-velocity  [m/s]  x-velocity  [m/s]  x-velocity  [m/s]  x-velocity  [m/s] 


time  =  100  s 


ure  4.3.  Comparison  of  analytical  and  computed  x-velocity  at  t 
(b)  200  s,  (c)  300  s,  and  (d)  400  s  for  r|0  =  0.025  m  in  Exai 
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4.2  Example  2  —  Salem  Harbor  Problem 

In  this  example,  we  used  the  model  to  simulate  hypothetical  contaminant  and  sediment 
transport  at  Salem  harbor  where  the  mean  sea  level  water  depth  varies  from  1 .4  m  to  8.4  m  [Yeh  and 
Kalasinsky,  1 977] .  We  first  used  the  flow  module  to  compute  a  quasi-steady  flow  pattern  controlled 
by  an  assumed  semidiurnal  tide.  Then  we  took  the  flow  result  at  the  end  of  a  tidal  cycle  as  the  initial 
condition  to  start  a  transient  simulation  where  a  dissolved  chemical  was  introduced  into  the  domain 
through  point  source  injection.  Twenty  semidiurnal  tidal  cycles  were  applied  to  perform  a  10  day 
simulation.  The  time  step  sizes  were  set  constant  at  20  seconds  for  the  flow  computation  and  60 
seconds  for  the  transport  computation. 

As  shown  in  Figure  4.4  for  the  flow  simulation,  we  had  a  tidal  boundary  condition  given  on 
the  ocean  boundary  side  with  a  maximum  tide  amplitude  of  1 .4  m  and  the  closed  boundary  condition 
applied  to  the  rest  of  the  boundary.  When  the  transient  simulation  began,  we  had  10  point  sources 
of  1  mVs  taken  into  account  to  simulate  man-induced  injection  (e.g.,  outfall).  Convergence  was 
considered  reached  when  either  the  maximum  relative  error  of  water  depth  was  less  than  10"4  or  the 
root  mean  square  error  was  less  than  10"8. 

The  transport  system  contained  (1)  3  suspended  sediments  (SSI,  SS2,  and  SS3)  and  3 
corresponding  bed  sediments  (BS1,  BS2,  and  BS3),  (2)  1  dissolved  chemical  species  in  the  water 
column  (C),  (3)  3  particulate  chemical  species  adsorbed  on  suspended  sediments  (PS1,  PS2,  and 
PS3;  one  for  each  sediment)  and  3  particulate  chemical  species  on  bed  sediments  (PB1,  PB2,  and 
PB3),  and  (4)  1  dissolved  chemical  species  in  the  interstitial  water  (CB).  The  following  adsorption 
reactions  were  taken  into  account. 


SSI 

kf  = 

0.001, 

SS2 

kf  = 

0.0005 

SS3 

kf  = 

0.0001 

BS1 

kf  = 

0.0001 

BS2 

kf  = 

0.0000 

BS3 

kf  = 

0.0000 

+  BS1 

kf  = 

0.0001 

+  BS2 

kf  = 

0.0000 

+  BS3 

kf  = 

0.0000 

where  kf  and  kb  represent  forward  and  backward  reaction  rate  constants,  respectively.  We  assumed 
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all  sediments  were  of  cohesive  type  and  deposition  and  erosion  rates  were  estimated  by  the  following 
equations  [Yeh,  1983]. 


S 


n 


if 


if 


-  1 


TcRn 


if' 

if  x 


where  Vsn  is  the  settling  velocity  of  the  n-th  size  fraction  sediment  [L/T] ;  TcDn  is  the  critical  shear 
stress  for  the  deposition  of  the  n-th  size  fraction  sediment  [M/L/T2];  TcRn  is  the  critical  shear  stress 
for  the  erosion  of  the  n-th  size  fraction  sediment  [M/L/T2];  En  is  the  erodibility  of  the  n-th  size 
fraction  sediment  [M/L2] ;  and  Tb  is  the  bottom  shear  stress  or  the  bottom  friction  stress  [M/L/T2]  that 
can  be  estimated  with  Manning’s  formula  [Chow,  1973].  Table  4.1  lists  the  characteristics  of  the 
three  sizes  of  sediments,  where  the  bulk  density  and  porosity  are  for  bed  sediments. 

Initially,  there  was  no  chemical  nor  bed  sediment  in  the  domain.  Only  the  three  suspended 
sediments  with  100  kg/m3  of  each  size  existed  throughout  the  region  of  interest.  The  water  volume 
exchange  rate  was  0.002  m3/m2/s  between  the  column  and  the  interstitial  waters.  The  flushing  decay 
rate  was  0.1  s"1  for  the  ocean  boundary  condition.  The  dispersion  coefficient  was  5.2  m2/s.  Neither 
radioactive  decay  nor  volatilization  existed.  Each  point  source  injected  the  dissolved  chemical  with 
a  rate  of  1  kg/s.  The  allowed  maximum  relative  error  of  concentration  was  set  to  10"4  for 
determining  convergence. 

Figure  4.5  plots  the  velocity  change  during  the  twentieth  tidal  cycle  of  the  transient 
simulation,  which  represents  a  quasi-steady  flow  pattern.  Based  on  the  computational  results,  this 
quasi-steady  flow  pattern,  determined  subject  to  both  tides  and  point  sources,  has  been  reached  after 
three  tidal  cycles.  Figures  4.6  and  4.7  plot  the  distributions  at  time  =  5  and  10  days  for  BS 1  and  SSI, 
respectively.  Figure  4.7  shows  the  trend  to  settle  SSI  in  the  bay  area  subject  to  the  periodic  tide. 
The  decrease  of  SSI  with  time  (Figure  4.7)  was  the  result  of  both  deposition  (Figure  4.6)  and  the 
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flushing  process  from  the  ocean  boundary  side.  Figures  4.8  and  4.9  show  the  distributions  of  BS2 
and  SS2,  respectively,  at  time  =  5  and  10  days.  The  concentration  contours  at  time  =  5  and  10  days 
are  barely  differentiable,  suggesting  that  a  quasi-steady  distribution  pattern  between  BS2  and  SS2 
may  be  reached  in  5  days.  An  examination  of  the  computational  results  revealed  that  a  quasi- steady 
distribution  pattern  can  be  considered  reached  after  5  tidal  cycles  (i.e.,  2.5  days).  Figure  4. 10  depicts 
the  contour  of  BS3  at  time  =  10  days,  which  presents  a  steady  state  concentration  distribution  that 
has  been  reached  in  this  10  day  simulation.  By  checking  the  numerical  results,  it  took  almost  no  time 
to  turn  all  SS3  into  BS3  due  to  its  large  settling  speed. 

Figure  4. 1 1  depicts  the  distributions  of  C  at  various  times.  It  is  obvious  that  the  periodic  flow 
pattern  could  effectively  transport  chemicals  from  the  injection  locations  to  the  northwestern  comer 
of  the  bay  through  advection.  Figure  4. 12  plots  the  contour  of  CB  at  time  =10  days,  which  matches 
that  of  C  (Figure  4. 1 1)  because  of  the  diffusive  exchange  of  chemicals  between  the  column  and  the 
interstitial  waters  as  well  as  the  identical  adsorption  reactions  onto  bed  sediments  for  C  and  CB  (the 
last  six  reactions  in  Eq.  (4.7)).  Figure  4.13  shows  the  distributions  of  PS1  and  PS2,  presented  as 
mass  of  particulate  chemicals  included  in  a  unit  volume  of  fluid  in  the  water  column  [kg/m3] ,  at  time 
=  10  days.  The  contour  patterns  are  also  consistent  with  that  of  C  because  the  distributions  of  SSI 
and  SS2  were  basically  smooth  over  the  bay  area  (Figures  4.7  and  4.9).  On  the  other  hand,  however, 
the  distributions  of  particulate  chemicals  on  bed  sediments,  presented  as  mass  of  particulate 
chemicals  per  unit  horizontal  bed  area  [kg/m2],  are  ragged  (Figure  4.14)  due  mainly  to  the  rough 
distributions  of  bed  sediments  (Figures  4.6,  4.8,  and  4.10). 
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Table  4.1.  List  of  sediment  parameters  in  example  2. 


SSI  and  BS1 

SS2  and  BS2 

SS3  and  BS3 

Settling  Speed  (m/s) 

1.2xl0‘6 

1.5xl0"4 

6.5xl0"2 

Critical  shear  stress  for  deposition  (kg/m/s2) 

0.5 

0.5 

0.6 

Critical  shear  stress  for  erosion  (kg/m/s2) 

0.4 

0.4 

0.5 

Erodibility  (kg/m2) 

1.0 

1.0 

1.0 

Specific  weight  (kg/m3) 

1307.7 

1500 

1727.3 

Bulk  density  (kg/m3) 

1200 

1300 

1400 

Porosity  (m3/m3) 

0.35 

0.4 

0.45 

Ocean 

Boundary 


Figure  4.4.  Discretization,  ocean  boundary,  and  point  sources  of  Salem  harbor 

in  the  demonstrative  example. 
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(a) 


(b) 


Figure  4.5.  Velocity  at  (a)  0.25T,  (b)  0.5T,  (c)  0.75T,  and  (d)  T  of  the  twentieth  tidal  cycle 
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(a) 

time  =  5  days 


Figure  4.6.  Numerical  results  of  BS1 


(a) 

time  =  5  days 


Figure  4.7.  Numerical  results  of  SSI 


(b) 

time  =  10  days 


[kg/m2]  at  time  =  (a)  5  days  and  (b)  10  days. 


(b) 

time  =  10  days 


[kg/m3]  at  time  =  (a)  5  days  and  (b)  10  days  . 
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(a)  (b) 

time  =  5  days  time  =  10  days 


Figure  4.8.  Numerical  results  of  BS2  [kg/m2]  at  time  =  (a)  5  days  and  (b)  10  days. 


Figure  4.9.  Numerical  results  of  SS2  [kg/m3]  at  time  =  (a)  5  days  and  (b)  10  days. 
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Figure  4.11.  Numerical  results  of  C  [kg/m3]  at  time  =  (a)  1  day,  (b)  3  days, 

(c)  5  days,  and  (d)  10  days. 


Figure  4.12.  Numerical  results  of  CB  [kg/m3]  at  time  =  10  days. 


Figure  4.13.  Numerical  results  of  (a)  PS1  [kg/m3]  and  (b)  PS2  [kg/m3]at  time  =  10  days. 
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(a) 

time  =  10  days 


(b) 


time  =  10  days 


Figure  4.14.  Numerical  results  of  (a)  PB1  [kg/m2],  (b)  PB2  [kg/m2], 
and  (c)  PB3  [kg/m2]  at  time  =  10  days. 
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4.3  Example  3  —  San  Diego  Bay  Problem  (I) 

In  this  example,  we  used  our  model  to  simulate  a  hypothetical  sediment  and  contaminant 
transport  at  San  Diego  South  bay  under  a  quasi-steady  flow  pattern  determined  according  to  an 
assumed  semidiurnal  tide.  This  quasi-steady  flow  pattern  was  computed  by  a  hydrodynamic  flow 
model  which  solved  shallow  water  equations  with  our  flow  module.  The  domain  was  discretized 
with  1415  elements  and  1567  nodes  (Figure  4.15).  The  mean  depth  varied  from  20  meters  near  the 
ocean  boundary  to  about  0.4  meter  toward  the  south  end  of  the  bay.  The  flow  pattern  was  determined 
with  a  tidal  boundary  condition  implemented  on  the  ocean  boundary  side  where  the  maximum  tidal 
amplitude  was  1.2  m  and  with  the  rest  of  the  boundary  treated  as  closed  (Figure  4.15).  It  was  also 
assumed  subject  to  8  point  sources  (Figure  4.15)  each  with  an  injection  rate  of  1  m3/s.  Manning’s 
n  was  assumed  0.01  throughout  the  entire  bay.  Figure  4.16  depicts  the  contours  of  water  surface 
elevation  in  the  bay  area  at  various  times  during  one  tidal  cycle. 

The  hypothetic  transport  system  contained  2  sizes  of  cohesive  sediments  and  1  dissolved 
chemical.  The  deposition  and  erosion  rate  were  computed  with  Eqs.  (4.8)  and  (4.9).  Table  4.2  lists 
the  characteristics  of  the  two  sizes  of  sediments,  where  the  bulk  density  and  porosity  were  used  for 
their  corresponding  bed  sediments.  Through  reactions,  we  can  have  a  total  of  4  sediments  (2 
suspended,  SSI  and  SS2,  and  2  bed,  BS1  and  BS2)  and  6  chemical  species  (one  dissolved  in  the 
water  column,  C,  one  dissolved  in  the  interstitial  water,  CB,  and  one  particulate  adsorbed  onto  each 
sediment,  PS1,  PS2,  PB1,  and  PB2)  included  in  our  transport  system.  The  following  adsorption 
reactions  were  considered. 

-  SSI  kf  =  0.0001 

-  SS2  kf  =  O.OOOf 

+  BS1  kf  =  0.0001 

+  BS2  kf  =  0.000^ 

+  BS1  kf  =  0.0001 
+  BS2  kf  =  0.0001 

Initially,  there  was  no  dissolved  chemical  nor  suspended  sediments  anywhere  in  the  domain. 
It  was  assumed  only  the  two  bed  sediments  of  100  kg/m2  were  uniformly  distributed  throughout  the 
region  of  interest.  As  the  simulation  began,  the  dissolved  chemical  was  released  into  the  domain 
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through  point  sources  with  a  rate  of  1  kg/s  at  each  point  source.  The  water  volume  exchange  rate  was 
0.002  m7m2/s  between  the  column  and  the  interstitial  waters.  The  flushing  decay  rate  was  0. 1  s"1  for 
the  ocean  boundary  condition.  The  dispersion  coefficient  was  5.2  m2/s.  Chemicals  were  neither 
radioactive  nor  volatile.  The  allowed  maximum  relative  error  of  concentration  was  10"4.  A  10  day 
simulation  was  performed  with  a  constant  time  step  size  of  20  seconds. 

Figures  4.17  and  4.18  show  the  distributions  of  BS1  and  BS2,  respectively,  at  time  =  5  and 
10  days.  It  is  observed  that  after  5  days  almost  all  bed  sediments  in  zone  A  were  washed  out,  while 
sediments  which  were  eroded  from  zone  A  were  pushed  into  zone  B  by  tides  and  then  deposited  due 
to  the  reduction  of  flow  velocities  along  the  bay’s  longitudinal  direction  from  the  ocean  boundary 
to  its  south  end.  In  zone  C  where  the  flow  velocity  was  too  small  to  resuspend  bed  sediments,  bed 
sediment  concentrations  did  not  change  much.  The  increase  of  bed  sediments  in  this  zone  was  due 
to  deposition  of  suspended  sediments  which  were  transported  by  tides  from  nearby  areas.  As  one 
can  see  from  Figure  4. 19  that  only  minor  amounts  of  suspended  sediments  were  transported  into  zone 
C. 

Figures  4.20  and  4.21  present  concentration  distributions  of  the  dissolved  chemical  in  the 
water  column  (i.e.,  C)  and  the  interstitial  water  (i.e.,  CB),  respectively,  at  both  time  =  5  days  and  time 
=  10  days.  Through  the  comparison  of  these  two  figures,  C  and  CB  had  approximately  the  same 
distribution  pattern.  In  addition,  their  distributions  did  not  change  drastically  from  time  =  5  days  to 
time  =10  days.  This  was  because  the  semidiurnal  tide  was  not  strong  enough  to  overcome  the  bay’ s 
long-strip  geometry  to  bring  out  chemicals  through  advection.  In  other  words,  dispersion  was  the 
main  process  to  bring  chemicals  away  from  the  injection  sources. 

Figures  4.22  and  4.23  plot  particulate  chemical  concentrations  on  suspended  sediments  and 
bed  sediments,  respectively,  at  time  =  10  days.  It  is  observed  that  there  existed  more  PB2  than  PB 1 
because  more  adsorption  was  allowed  for  dissolved  chemicals  to  be  adsorbed  onto  BS2  by  comparing 
the  last  four  reactions  inEq.  (4.10).  On  the  other  hand,  we  observed  more  PS  1  thanPS2  even  though 
the  system  should  allow  more  PS2  to  be  formed  according  to  the  first  two  adsorption  reactions  in  Eq. 
(4.10).  This  was  because  SSI  was  at  least  one  order  of  magnitude  larger  than  SS2  and  we  present 
PS  1  and  PS2  in  Figure  4.23  as  the  mass  of  particulate  chemicals  on  bed  sediments  per  unit  bed  area. 
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Table  4.2.  List  of  sediment  parameters  in  Example  3. 


Sediment  1 
(SSI  and  BS1) 

Sediment  2 
(SS2  and  BS2) 

Settling  Speed  (m/s) 

1.2xl0"6 

1.5xl0"4 

Critical  shear  stress  for  deposition  (kg/m/s2) 

0.5 

0.5 

Critical  shear  stress  for  erosion  (kg/m/s2) 

0.4 

0.4 

Erodibility  (kg/m2) 

1.0 

1.0 

Specific  weight  (kg/m3) 

1500 

1285.7 

Bulk  density  (k/m3) 

1300 

1200 

Porosity  (m3/m3) 

0.4 

0.3 
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Figure  4.16.  Surface  elevation  at  (a)  0  or  T,  (b)  0.25T,  (c)  0.5T,  and  (d)  0.75T 
in  the  bay  area  during  a  tidal  cycle  where  T  =  12  hours. 
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Figure  4.17.  Numerical  results  of  BS1  [kg/m2]  at  time  =  (a)  5  days  and  (b)  10  days 

in  the  bay  of  Example  4. 


Figure  4.18.  Numerical  results  of  BS2  [kg/m2]  at  time  =  (a)  5  days  and  (b)  10  days 

in  the  bay  of  Example  3. 
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Figure  4.19.  Numerical  results  of  (a)  SSI  [kg/m3]  and  (b)  SS2  [kg/m3] 
at  time  =10  days  in  the  bay  of  Example  3. 
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Figure  4.20.  Numerical  results  of  C  [kg/m3]  at  time  =  (a)  5  days  and  (b)  10  days 

in  the  bay  of  Example  3. 


Figure  4.21.  Numerical  results  of  CB  [kg/m3]  at  time  =  (a)  5  days  and  (b)  10  days 

in  the  bay  of  Example  3. 
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Figure  4.22.  Numerical  results  of  (a)  PB1  [kg/m2]  and  (b)  PB2  [kg/m2] 
at  time  =10  days  in  the  bay  of  Example  3. 


Figure  4.23.  Numerical  results  of  (a)  PS1  [kg/m3]  and  (b)  PS2  [kg/m3] 
at  time  =  10  days  in  the  bay  of  Example  3. 
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4.4  Example  4  —  San  Diego  Bay  Problem  (II) 

This  example  differed  from  Example  3  only  in  the  size  of  the  simulated  region  and  the 
application  of  open  boundary  conditions.  As  shown  in  Figure  4.24  (a),  the  domain  was  discretized 
with  2467  elements  and  2566  nodes  where  the  grids  used  to  discretize  the  bay  area  were  identical  to 
those  used  in  Example  3.  The  mean  depth  varied  from  about  85  meters  close  to  the  Point  C  on  the 
ocean  boundary  to  about  0.4  meter  toward  the  south  end  of  the  bay  (Figure  4.24(b)).  A  tidal 
boundary  condition  with  the  maximum  tidal  amplitude  of  2.0  m  was  applied  to  the  ocean  boundary 
where  a  wave  lead  of  4  s  per  every  100  m  was  assumed  from  Point  A  to  B,  A  to  C,  and  C  to  D,  i.e., 
the  wave  came  from  northeast  and  arrived  at  Point  A  first,  then  Point  B,  Point  C,  and  finally  Point 
D.  The  rest  of  the  boundary  was  treated  as  closed  as  that  in  Example  3(Figure  4.24). 

Figure  4.25  depicts  the  contours  of  water  surface  elevation  at  various  times  during  one  tidal 
cycle.  It  is  observed  that  Figures  4.25  and  4.16  show  similar  elevation  distributions  in  the  bay  at 
various  times.  Therefore,  we  also  expected  similar  concentration  distributions  of  the  dissolved 
chemical  in  the  water  column  (i.e.,  C)  and  the  interstitial  water  (i.e.,  CB)  in  both  Examples  3  and  4 
because  the  sources  were  injected  in  the  middle  of  the  bay  and  the  tide  could  not  transport  the 
chemical  to  the  Pacific  Ocean  through  the  exit  of  the  bay  in  10  days.  Figures  4.26  and  4.27  present 
concentration  distributions  of  C  and  CB,  respectively,  at  both  time  =  5  days  and  time  =10  days. 
Through  the  comparison  of  these  two  figures  with  Figures  4.20  and  4.2 1 ,  we  verified  our  expectation. 

Figures  4.28  and  4.29  show  the  distributions  of  BS1  and  BS2,  respectively,  at  time  =  5  and 
10  days.  Fike  the  chemical  concentrations  mentioned  above,  similar  numerical  results  were  obtained 
in  the  bay,  except  for  the  throat  part  that  suspended  sediments  may  enter  the  bay  from  or  exit  to  the 
Pacific  Ocean.  Figure  4.30  helps  to  demonstrate  this  by  showing  the  distribution  of  suspended 
sediments  at  time  =10  days  where  suspended  sediments  in  the  throat  area  were  more  than  those  in 
Example  3  (Figure  4.19). 

Figures  4.31  and  4.32  plot  the  distribution  of  particulate  chemicals  at  time  =10  days.  Again, 
the  distributions  similar  to  those  in  Example  3  (Figures  4.22  and  4.23).  It  is  thus  suggested  a 
simulated  domain  that  includes  only  the  bay  can  be  used  to  save  computer  time  when  the  injection 
is  placed  far  enough  from  the  throat  of  the  bay. 
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Figure  4.25.  Surface  elevation  at  (a)  0  or  T,  (b)  0.25T,  (c)  0.5T,  and  (d)  0.75T 
during  a  tidal  cycle  in  Example  4  where  T  =  12  hours. 


r- 

(N 


Figure  4.27.  Numerical  results  of  CB  [kg/m3]  at  time  =  (a)  5  days  and  (b)  10  days  in  Example  4. 
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Figure  4.32.  Numerical  results  of  (a)  PS1  [kg/m3]  and  (b)  PS2  [kg/m3]  at  time  =  10  days  in  Example  4. 
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Chapter  5.  SUMMARY 

We  have  presented  a  2-D  depth-averaged  model  simulating  hydrodynamics  and  reactive 
contaminant  and  sediment  transport  in  bay /estuary  areas.  A  splitting  strategy  was  employed  in  the 
flow  module  to  solve  shallow  water  equations.  We  used  the  method  of  characteristics  first  to  solve 
the  equations  without  including  eddy  flux  terms  where  local  diagonalization  was  achieved,  and  then 
we  accounted  for  these  terms  with  the  Galerkin  finite  element  method.  Zero  normal  flux  was  applied 
to  the  closed  boundary,  while  zero  lateral  flux  was  assumed  on  the  open  boundary.  Both  were  to 
replace  the  first  characteristic  equations  of  boundary  nodes.  Either  normal  flux-  or  water  surface- 
driven  boundary  conditions  can  be  used  on  the  open  boundary  where  sub-critical  flow  exists.  The 
open  boundary  condition  was  to  replace  the  equation  of  the  incoming  characteristic,  either  the  second 
or  the  third,  of  boundary  nodes.  We  employed  a  generic,  mechanistic  approach  to  model  reactive 
contaminant  and  sediment  migration  in  the  transport  module,  desiring  to  promote  reactive 
contaminant  transport  modeling  in  this  field  so  that  practical  systems  can  be  more  accurately 
simulated.  Chemical  reactions,  including  aqueous  complexation,  adsorption/desorption,  and 
volatilization,  were  assumed  elementary  and  kinetic  where  reaction  rates  were  computed  through 
forward  and  backward  reaction  rate  constants  based  on  the  collision  theory.  Suspended  sediments, 
dissolved  chemicals  in  the  water  column,  and  particulate  chemicals  on  suspended  sediments  were 
subject  to  transport  and  thus  were  mobile.  Bed  sediments,  dissolved  chemicals  in  the  interstitial 
water,  and  particulate  chemicals  on  bed  sediments  were  considered  immobile.  The  set  of  governing 
equations  were  solved  with  a  predictor-corrector  strategy.  In  the  predictor  step,  the  physical  transport 
terms  were  solved  by  using  the  Lagrangian-Eulerian  finite  element  method  with  all  source/sink  terms 
evaluated  at  the  previous  time.  In  the  corrector  step  where  non-linearity  was  handled,  the  Picard 
method  was  used  to  solve  for  sediment  distributions,  and  the  Newton-Raphson  method  was 
employed  to  solve  the  correction  equations  of  reactive  contaminants.  An  adaptive  explicit-implicit 
scheme  was  also  applied  to  avoid  negative  numerical  results.  Four  example  problems  accounting 
for  1-D  standing  wave  problem  and  the  transports  of  both  sediments  and  contaminants  in  the  Salem 
harbor  as  well  as  in  the  Sandiego  South  bay  with  semidiurnal  tides  applied  on  the  ocean  boundary 
were  used  to  demonstrate  the  capability  of  the  model. 
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Appendix  I.  DATA  INPUT  GUIDE 

This  data  input  guide  describe  the  following  files:  geometry  file,  flow  model  file, 
chemical/sediment  and  reaction  information  file,  transport  model  file,  initial  water  surface  elevation 
file,  initial  velocity  file,  initial  dissolved  chemical  concentration  file,  initial  particulate  chemical 
concentration  file,  and  initial  sediment  concentration  file. 


1.1.  File  1:  geometry  file  (GEOM  file) 

T1-T3  Cards  Job  title 

Only  one  Tl,  T2,  and  T3  card  can  be  used. 

GE  Card  Element  information 


Field 

Variable 

Value 

Description 

C  1-2 

IC1 

GE 

Card  group  identifier. 

C  3 

IC3 

3 

Identification  of  element. 

A  triangular  element  is  being  input. 

4 

A  quadrilateral  element  is  being  input. 

1 

M 

+ 

id  of  global  element. 

2 

IEl(M.l) 

+ 

id  of  the  global  node  that  serves  as  the  1-st  element  node  of  the 
global  element. 

3 

IE2(M,2) 

+ 

id  of  the  global  node  that  serves  as  the  2-nd  element  node  of  the 
global  element.. 

4 

IE2(M,3) 

+ 

id  of  the  global  node  that  serves  as  the  3-rd  element  node  of  the 
global  element. 

5 

IE2(M,4) 

+ 

id  of  the  global  node  that  serves  as  the  4-th  element  node  of  the 
global  element.  It  is  zero  for  a  triangular  element. 

6 

IE2(M,5) 

+ 

id  of  the  material  type  of  the  global  element. 

GN  Card 

Node  information 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

GN 

Card  group  identifier. 

C  3 

IC3 

blank 

1 

N 

+ 

id  of  global  node.. 

2 

X2(N,1) 

+ 

The  x-coordinate  of  the  global  node. 

3 

X2(N,2) 

+ 

The  y-coordinate  of  the  global  node. 

4 

X2(N,3) 

+ 

The  mean  water  surface  elevation  of  the  global  node. 

EN  Card 

End  of  data  control 

Field 

Variable 

Value 

Description 

Cl -2 

IC1 

EN 

Card  group  identifier. 

C  3 

IC3 

D 

End  of  input  data. 
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1.2.  File  2:  flow  file  (BCFT  file) 


T1-T3  Cards  Job  title 

Only  one  Tl,  T2,  and  T3  card  can  be  used. 


OP  Card 
Field  Variable 
C  1-2  IC1 
C  3  IC3 
1  KMOD 


2  IRUN 


Run  option  parameters 


Value 

OP 

1 


0 

1 

2 

12 


Description 

Card  group  identifier. 

Simulation  selection. 

Index  of  flow/transport  simulation. 

Neither  flow  nor  transport  simulation  will  proceed.  Only  data  are 
read  for  testing  purpose. 

Only  flow  is  simulated. 

Only  transport  is  simulated. 

Both  flow  and  transport  are  simulated. 

Index  of  flow  simulation  run. 

Quasi-steady  flow  is  input  from  the  file  “STEADY.STO”  for 
executing  transient  transport  simulations. 

Initial  flow  conditions  are  given  for  transient  simulations. 
Quasi-steady  flow  will  be  compute  under  given  boundary 
conditions. 

Quasi-steady  flow  will  be  computed  first  and  then  be  used  to 
simulate  transient  transport. 

Quasi-steady  flow  will  be  computed  first  and  the  final  result  will 
be  used  as  the  initial  condition  for  transient  simulations. 


KMOD 

IRUN 

Simulation 

0 

any 

Only  input  data  are  read 

1 

-1 

Transient  flow 

1 

0 

Quasi-steady  flow 

1 

2 

Quasi-steady  flow  &  transient  flow 

2 

-2 

Transient  transport 

12  or  2 

1 

Quasi-steady  flow  and  transient  transport 

12 

-1 

Transient  flow  and  transport 

12 

2 

Quasi-steady  flow  &  transient  flow  and  transport 

*  (1,0)  +  (2,-2)  =  (12,1)  or  (2,1) 
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Field 

Variable 

Value 

Description 

C  1-2 

IC1 

OP 

Card  group  identifier. 

C  3 

IC3 

2 

Solver  specification. 

1 

ILUMP 

0 

Index  of  using  mass  lumping. 

No  mass  lumping. 

1 

Mass  lumping. 

2 

IPNTS 

1 

Index  of  using  linear  matrix  solver  for  flow  simulations. 

The  pointwise  iterative  matrix  solver. 

2 

Preconditioned  Conjugate  Gradient  Method  (polynomial), 

3 

Preconditioned  Conjugate  Gradient  Method  (Incomplete 
Cholesky). 

3 

IQUAR 

11 

Index  of  using  quadrature  for  numerical  integration. 

Nodal  quadrature  for  both  line  and  element  integrations. 

12 

Nodal  quadarature  for  line  integration  and  Gaussian  quadrature  for 
element  integration. 

21 

Gaussian  quadarature  for  line  integration  and  nodal  quadrature  for 
element  integration. 

22 

Gaussian  quadrature  for  both  line  and  element  integrations. 

4 

IDIFF 

1 

Index  of  taking  into  account  the  eddy  viscosity  terms. 

Yes. 

0 

No,  the  eddy  viscosity  terms  are  neglected. 

5 

[ALLOW 

1 

Index  of  allowing  the  simulation  to  continue  when  there  is  no 
convergence  reached  for  the  non-linear  iteration. 

Yes. 

0 

No. 

6 

1MODEL 

1 

Index  of  simulating  flow. 

Using  MOC  to  solve  the  dynamic  model  only. 

-10 

Specifying  water  depth  and  flow  velocity  with  the  initial  condition. 
The  flow  equation  is  not  solved. 

7 

IQUASI 

0 

Index  of  quasi-steady  state  check. 

No  check. 

1 

Performing  the  quasi-steady  check 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

OP 

Card  group  identifier. 

C  3 

IC3 

3 

Weighting  factor  for  simulations. 

1 

W 

0.5 

Time  derivative  weighting  factor  for  flow  simulation. 
Crank-Nicolson  central. 

1.0 

Backward  difference. 

2 

OME 

0.0- 1.0 

Relaxation  parameter  for  solving  nonlinear  matrix  equations. 
Under  relaxation. 

1.0 

Exact  relaxation. 

1. 0-2.0 

Over  relaxation. 

3 

OM1 

0.0- 1.0 

Relaxation  parameter  for  solving  linearized  matrix  equation. 
Under  relaxation. 

1.0 

Exact  relaxation. 

1. 0-2.0 

Over  relaxation. 

1-4 


Field 

Variable 

Value 

Description 

C  1-2 

IC1 

OP 

Card  group  identifier. 

C  3 

IC3 

4 

Preconditioned  Conjugate  Gradient  method. 

1 

IEIGEN 

0 

GG  will  not  be  read. 

1 

GG  will  be  read. 

2 

GG 

+ 

Upper  bound  of  the  maximum  eigenvalue  of  the  coefficient  matrix 
used  in  the  Preconditioned  Conjugate  Gradient  method. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

OP 

Card  group  identifier. 

C  3 

IC3 

5 

Method  of  characteristics. 

1 

ISK 

-1 

Index  of  choosing  the  directions  of  characteristics. 

Using  the  velocity  direction  as  the  direction  of  characteristics. 

0 

Solve  the  diagonalization  equations  for  the  directions. 

1 

Specify  directions  by  users. 

2 

ANGLE  1 

+ 

User  specified  angle  for  the  direction  of  the  first  characteristic. 
(COS(ANGLEl),  SIN(ANGLEl))  is  the  direction. 

3 

ANGLE2 

+ 

User  specified  direction  for  the  second  and  third  characteristics. 
(COS(ANGLE2),  SIN(ANGLE2))  is  the  direction. 

IP  Card 

Iteration  Parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IP 

Card  group  identifier. 

C  3 

IC3 

1 

Total  simulation  time. 

1 

NNITER 

+ 

Allowed  number  for  the  non-linear  iteration  in  the  flow 
simulation. 

2 

NPITER 

+ 

Allowed  number  for  the  linearized  iteration  in  the  flow 
simulation. 

3 

TOLBH 

+ 

Allowed  relative  error  for  water  depth. 

4 

TOLBV 

+ 

Allowed  relative  error  for  velocity  components. 

5 

V_CUT 

+ 

Cut-off  of  velocity.  If  the  computed  velocity  is  smaller  than  this 
value,  it  is  set  to  zero. 

6 

FMSR 

+ 

Multiplication  factor  to  set  the  allowed  mean  square  root  errors, 
BHMSR  for  water  depth  and  B  VMSR  for  velocity,  where  BHMSR 
=  TOLBH*FMSR;  BVMSR  =  TOLBV*FMSR. 

7 

ACCU 

+ 

The  allowed  absolute  error  between  two  tidal  cycles  to  reach  quasi¬ 
steady  state. 

PT  Card 

Time  inarching  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

PT 

Card  group  identifier. 

C  3 

IC3 

1 

Total  simulation  time. 

1 

NXW 

+ 

Subelement  number  in  one  direction  that  is  used  in  the  “in- 
element”  particle  tracking.  The  total  subelement  number  included 
in  one  global  element  is  NXW*NXW. 

2 

IDETQ 

1 

Index  of  particle  tracking  approach. 

Using  the  average-velocity  approach. 

2 

Using  the  constant-velocity  approach. 
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TC  Card 

Iteration  Parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

TC 

Card  group  identifier. 

C  3 

IC3 

1 

Total  simulation  time. 

1 

TMAX 

+ 

Maximum  simulation  time.  [T] 

2 

DELT 

+ 

Time  duration  (period)  of  one  tidal  cycle.  [T] 

3 

NTI 

+ 

No.  of  tidal  cycles. 

4 

NTIC 

+ 

No.  of  time  steps  in  a  tidal  cycle. 

5 

NTCK 

+ 

No.  of  time  step  intervals  to  check  quasi-steady  flow.  For  instance, 
if  NTCK  =  10,  we  check  quasi-steady  flow  by  comparing  the 
computational  results  every  10  time  steps  from  corresponding 
times  in  successive  tidal  cycles. 

6 

NTSTO 

+ 

No.  of  time  step  intervals  to  record  flow  results  for  future  use  after 
the  quasi-steady  condition  is  reached. 

7 

NTRANSP 

+ 

No.  of  time  step  intervals  to  compute  transport.  For  instance,  if 
NTRANSP  =  5,  we  perform  one  transport  computation  in  five  flow 
time  steps. 

OC  Card 

Output  control  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

OC 

Card  group  identifier. 

C  3 

IC3 

1 

Print  control  parameters. 

1 

JOPT 

0 

Print  at  specified  intervals. 

1 

Print  at  specified  set  of  time  values. 

2 

KPRT 

+ 

Specify  the  interval  if  JOPT  =  0. 

NPRINT 

+ 

Specify  total  number  of  the  specified  times,  if  JOPT=l . 

Note: 

NPRINT  new  lines  after  OC1  card  are  needed  if  JOPT  =  1. 

(the  first  new  line) 

1 

TPRT(l) 

+ 

The  1-st  specified  time  moment  at  which  the  associated  simulation 
results  are  to  be  printed. 

(the  second  new  line) 

1 

TPRT(2) 

+ 

The  2-nd  specified  time  moment  at  which  the  associated  simulation 
results  are  to  be  printed. 

(the  NPRINT  new  line) 

1 

TPRT(NPRINT) 

+ 

The  NPRINT -th  specified  time  moment  at  which  the  associated 
simulation  results  are  to  be  printed. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

OC 

Card  group  identifier 

C  3 

IC3 

2 

Print  options. 

1 

NSELT 

+ 

Total  options  to  be  selected. 

2 

KPRO(i) 

1 

Print  options. 

Print  water  surface  elevation. 

2 

Print  velocity. 

3 

Print  water  depth. 
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Field 

Variable 

Value 

Description 

C  1-2 

IC1 

OC 

Card  group  identifier. 

C  3 

IC3 

3 

Output  control  for  post-process. 

1 

IFILE 

0 

Save  as  ascii  format. 

1 

Save  as  binary  format. 

2 

KOPT 

0 

Save  at  specified  interval. 

1 

Save  at  specified  set  of  times. 

3 

KDSK 

+ 

Save  every  specified  interval,  if  KOPT=0. 

NPOST 

+ 

Save  at  specified  set  of  times,  if  KOPT=l . 

Note: 

NPOST  new  lines  after  OC3  card  are  needed  if  KOPT  =  1 . 

(the  first  new  line) 

TSV(l) 

+ 

The  1  -st  specified  time  moment  at  which  the  associated  simulation 
results  are  to  be  saved. 

(the  second  new  line) 

TSV(2) 

+ 

The  2-nd  specified  time  moment  at  which  the  associated  simulation 
results  are  to  be  saved. 

(the  NPOST  new  line) 

TSV(NPOST) 

+ 

The  NPOST-th  specified  time  moment  at  which  the  associated 
simulation  results  are  to  be  saved. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

OC 

Card  group  identifier. 

C  3 

IC3 

4 

Solution  file  output  options. 

1 

KSELT 

+ 

Total  options  to  be  selected. 

2 

KSAVE(i) 

1 

Save  options. 

Save  water  surface  elevation. 

2 

Save  flux. 

3 

Save  water  depth. 

MP  Card 

Material 

property  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

MP 

Card  group  identifier. 

C  3 

IC3 

1 

Material  property. 

1 

i 

+ 

id  of  material  type. 

2 

PROP(l.i) 

+ 

Manning’s  roughness  associated  with  the  material  type. 

3 

PROP(2,i) 

+ 

xx-eddy  viscosity  associated  with  the  material  type.  [L2/T] 

4 

PROP(3,i) 

+ 

yy-eddy  viscosity  associated  with  the  material  type.  [L2/T] 

5 

PROP(4,i) 

+ 

xy-eddy  viscosity  associated  with  the  material  type.  [L2/T] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

MP 

Card  group  identifier. 

C  3 

IC3 

2 

Density  of  water  and  acceleration  of  gravity. 

1 

ITUNIT 

Index  of  time  unit. 

1  The  time  unit  for  simulation  is  second. 

2  The  time  unit  for  simulation  is  minute. 
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3 

The  time  unit  for  simulation  is  hour. 

2 

ILUNIT 

1 

Index  of  length  unit. 

The  length  unit  for  simulation  is  meter. 

2 

The  length  unit  for  simulation  is  foot. 

3 

The  length  unit  for  simulation  is  kilometer. 

4 

The  length  unit  for  simulation  is  mile. 

3 

RHO 

+ 

Density  of  water,  (M/L3). 

RF  Card 

Rainfall  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

RF 

Card  group  identifier. 

C  3 

IC3 

1 

Rainfall  data  for  the  simulation. 

1 

i 

+ 

id  of  global  element. 

2 

IRTYP(i) 

+ 

id  of  xy  series  to  describe  the  time -dependent  rainfall  rate  for  the 
element. 

SS  Card 

Source/sink  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

SS 

Card  group  identifier. 

C  3 

IC3 

1 

Source/sink  data  for  the  simulation. 

1 

i 

+ 

id  of  the  global  node. 

2 

ISR(i) 

+ 

id  of  xy  series  to  describe  the  time-dependent  source/sink  rate  for 

the  node. 


OB  Card 

Open  boundary  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

OB 

Card  group  identifier. 

C  3 

IC3 

1 

Dirichlet  boundary  nodes  for  the  simulation. 

1 

NPDB(j) 

+ 

id  of  global  node  corresponding  to  the  j-th  Dirichlet  boundary 
node. 

2 

IDTYP(j) 

+ 

id  of  xy  series  to  describe  the  time-dependent  water  surface 
elevation  or  water  depth  profile  associated  with  the  Dirichlet 
boundary  node. 

3 

IDBH(j) 

1 

Index  of  the  profile  type  for  the  j-th  Dirichlet  boundary  node. 
This  profile  is  a  time-dependent  water  surface  elevation  profile. 

-1 

This  profile  is  an  analytical  function  as  the  form  of  asin27I(t+t0)/T. 
The  associated  xy  series  needs  two  points  to  describe  a  and  t0. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

OB 

Card  group  identifier. 

C  3 

IC3 

2 

Flux-type  boundary  sides  for  the  simulation. 

1 

NSCB(j,l) 

+ 

id  of  global  node  corresponding  to  the  1-st  node  of  the  j-th  flux- 
type  boundary  side. 

2 

NSCB(j,2) 

+ 

id  of  global  node  corresponding  to  the  2-nd  node  of  the  j-th  flux- 
type  boundary  side. 

3 

ICTYP(j,l) 

+ 

id  of  xy  series  to  describe  the  time-dependent  normal  flux  profile 

associated  with  the  open  boundary  side. 
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IC  Card 

Initial  condition  data  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IC 

Card  group  identifier. 

C  3 

IC3 

H 

Water  surface  elevation. 

1 

IHEAD 

0 

Constant  water  surface  elevation. 

1 

Variable  water  surface  elevation.  In  this  case,  the  user  needs  to 
specify  an  ICWD  file  in  the  super  file  as  explained  in  the  end  of 
this  appendix. 

2 

HCONST 

+ 

Value  of  the  constant  water  surface  elevation  for  IHEAD  =  0.  [L] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IC 

Card  group  identifier. 

C  3 

IC3 

U 

Velocity. 

1 

IV 

0 

Constant  velocity. 

1 

Variable  velocity.  In  this  case,  the  user  needs  to  specify  an  ICQL 
file  in  the  super  file  as  explained  in  the  end  of  this  appendix. 

2 

UCONST 

+ 

Value  of  the  constant  velocity.  [L/T] 

3 

VCONST 

+ 

Value  of  the  constant  velocity.  [L/T] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IC 

Card  group  identifier. 

C  3 

IC3 

s 

Initial  condition  start  type. 

1 

ISTART 

0 

Cold  start,  water  depth/water  surface  elevation  and  velocity. 

1 

Hot  start,  water  surface  elevation  and  velocity. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IC 

Card  group  identifier. 

C  3 

IC3 

T 

Initial  condition  start  time. 

1 

HSTIME 

+ 

Time  corresponding  for  the  hot  start. 

XY  Card 

X-Y  Series  Parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

XY 

Card  group  identifier. 

C  3 

IC3 

1 

Generation  of  any  x-y  series  function. 

1 

i 

+ 

Index  number  for  x-y  series. 

2 

NPOINT(i) 

+ 

The  number  of  x-y  values  in  the  series. 

3 

DELTAX 

0 

Dummy  values. 

4 

DELTAY 

0 

Dummy  values. 

5 

REPEAT 

0 

Dummy  values. 

6 

BEGCYC 

0 

Dummy  values. 

7 

TNAME 

+ 

A  character  string  representing  the  name  of  the  XY  series. 

(new  line)X,Y 

After  XY  card,  the  x-y  values  of  the  series  are  listed  one  pair  per  line  up 

to  NPOINT(i). 

EN  Card 

End  of  data  control 

Field  Variable 

Value 

Description 

Cl -2  IC1 

EN 

Card  group 
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C  3  IC3  D  End  of  input  data. 


1.3.  File  3:  transport  file  (BCTT  file) 

T1-T3  Cards  Job  title 

Only  one  Tl,  T2,  and  T3  card  can  be  used. 

OP  Card 
Field  Variable 
C  1-2  IC1 
C  3  IC3 
1  IDCHEM 


2  IDSED 


Value  Description 

OP  Card  group  identifier. 

2  Solver  specification. 

Index  of  using  mass  lumping  for  transport  simulations. 

0  No  mass  lumping. 

1  Mass  lumping. 

Index  of  using  linear  matrix  solver  for  transport  simulations. 

1  The  pointwise  iterative  matrix  solver. 

2  Preconditioned  Conjugate  Gradient  Method  (polynomial), 

3  Preconditioned  Conjugate  Gradient  Method  (Incomplete 
Cholesky). 

Index  of  using  quadrature  for  numerical  integration. 

1 1  Nodal  quadrature  for  both  line  and  element  integrations. 

1 2  Nodal  quadarature  for  line  integration  and  Gaussian  quadrature  for 
element  integration. 

2 1  Gaussian  quadarature  for  line  integration  and  nodal  quadrature  for 
element  integration. 

22  Gaussian  quadrature  for  both  line  and  element  integrations. 

Index  of  using  adapted  implicit-explicit  scheme  for  the  predictor 

step. 

0  No. 

1  Yes. 


Field  Variable 
C  1-2  IC1 
C  3  IC3 
1  WT 


Value  Description 

OP  Card  group  identifier. 

3  Weighting  factor  for  transport  simulations. 

Time  derivative  weighting  factor  for  transport  simulation. 
0.5  Crank-Nicolson  central. 

1.0  Backward  difference. 


Field  Variable 
C  1-2  IC1 
C  3  IC3 
1  ILUMPT 


2  IPNTST 


3  IQUART 


4  IADAPT 


tsuii  option  paiameieis 


Value  Description 

OP  Card  group  identifier. 

1  Simulation  selection. 

Index  of  chemical  transport  simulation. 
0  Chemical  transport  is  not  simulated. 

1  Chemical  transport  is  simulated. 

Index  of  sediment  transport  simulation. 
0  Sediment  transport  is  not  simulated. 

1  Sediment  transport  is  simulated. 
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2 

OMET 

0.0- 1.0 

Relaxation  parameter  for  solving  nonlinear  matrix  equations. 
Under  relaxation. 

1.0 

Exact  relaxation. 

1. 0-2.0 

Over  relaxation. 

3 

OMIT 

0.0- 1.0 

Relaxation  parameter  for  solving  linearized  matrix  equation. 
Under  relaxation. 

1.0 

Exact  relaxation. 

1. 0-2.0 

Over  relaxation. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

OP 

Card  group  identifier. 

C  3 

IC3 

4 

Preconditioned  Conjugate  Gradient  method. 

1 

IEIGENT 

0 

GGT  will  not  be  read. 

1 

GGT  will  be  read. 

2 

GGT 

+ 

Upper  bound  of  the  maximum  eigenvalue  of  the  coefficient  matrix 
used  in  the  Preconditioned  Conjugate  Gradient  method. 

IP  Card 

Iteration  Parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IP 

Card  group  identifier. 

C  3 

IC3 

1 

Total  simulation  time. 

1 

NNITERT 

+ 

Allowed  number  for  the  non-linear  iteration. 

2 

NPITER 

+ 

Allowed  number  for  the  linearized  iteration. 

3 

TOLBC 

+ 

Allowed  relative  error  for  getting  convergent  transport  solutions 
for  chemicals. 

4 

TOLBS 

+ 

Allowed  relative  error  for  getting  convergent  transport  solutions 
for  sediments. 

PT  Card 

Particle  tracking  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

PT 

Card  group  identifier. 

C  3 

IC3 

1 

Particle  tracking  controlling  parameters. 

1 

NXWT 

+ 

Subelement  number  in  one  direction  that  is  used  in  the  “in¬ 
element”  particle  tracking.  The  total  subelement  number  included 
in  one  global  element  is  NXWT*NXWT. 

2 

IDETQT 

1 

Index  of  particle  tracking  approach. 

Using  the  average-velocity  approach. 

2 

Using  the  constant-velocity  approach. 

MP  Card 

Material  property  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

MP 

Card  group  identifier. 

C  3 

IC3 

1 

Material  property  for  transport. 

1 

i 

+ 

id  of  material  type. 

2 

PROPT(l,i) 

+ 

Longitudinal  dispersivity  associated  with  the  material  type.  [L] 

3 

PROPT(2,i) 

+ 

Lateral  dispersivity  associated  with  the  material  type.  [L] 

4 

PROPT(3,i) 

+ 

Pure  diffusion  coefficient  associated  with  the  material  type.  [L2/T] 

Ml 


SS  card 

Source/sink  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

SS 

Card  group  identifier. 

C  3 

IC3 

1 

Source/sink  profile  for  dissolved  chemicals. 

1 

NP 

+ 

id  of  the  global  node. 

2 

j 

+ 

id  of  dissolved  chemical. 

3 

ITSR(NP,j) 

+ 

id  of  xy  series  to  describe  the  time-dependent  source/sink  profile 
for  the  dissolved  chemical  at  the  global  node. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

SS 

Card  group  identifier. 

C  3 

IC3 

2 

Source/sink  profile  for  suspended  sediment. 

1 

NP 

+ 

id  of  the  global  node. 

2 

j 

+ 

id  of  suspended  sediment  size. 

3 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

4 

ITSR(NP.m) 

+ 

id  of  xy  series  to  describe  the  time-dependent  source/sink  profile 
for  the  suspended  sediment  size  at  the  global  node,  where  m  = 
NCHEM  +  j. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

SS 

Card  group  identifier. 

C  3 

IC3 

3 

Source/sink  profile  for  particulate  chemicals  on  suspended 
sediments. 

1 

NP 

+ 

id  of  the  global  node. 

2 

j 

+ 

id  of  suspended  sediment. 

3 

k 

+ 

id  of  particulate  chemical. 

4 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

5 

NSSIZE 

+ 

Number  of  sediment  sizes. 

6 

ITSR(NP.m) 

+ 

id  of  xy  series  to  describe  the  time-dependent  source/sink  profile 
for  the  particulate  chemical  on  the  suspended  sediment  size  at  the 
global  node,  where  m  =  NCHEM  +  NSSIZE  +  (j-1  )*NCHEM  +  k, 
j  E  [1, NSSIZE],  and  k  E[l,  NCHEM]. 

DB  Card 

Dirichlet 

boundary  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

DB 

Card  group  identifier. 

C  3 

IC3 

1 

Dirichlet  boundary  conditions  for  dissolved  chemicals. 

1 

NP 

+ 

id  of  Dirichlet  boundary  node. 

2 

NPDBT(NP) 

+ 

id  of  global  node  corresponding  to  the  Dirichlet  boundary  node. 

3 

i 

+ 

id  of  dissolved  Chemical. 

4 

IDTYPT(NP,i) 

+ 

id  of  xy  series  to  describe  the  of  time-dependent  concentration 
profile  associated  with  the  Dirichlet  boundary  node  for  the 
dissolved  chemical. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

DB 

Card  group  identifier. 

C  3 

IC3 

2 

Dirichlet  boundary  conditions  for  suspended  sediments. 

1 

NP 

+ 

id  of  Dirichlet  boundary  node. 
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2 

NPDBT(NP) 

+ 

id  of  global  node  corresponding  to  the  Dirichlet  boundary  node. 

3 

k 

+ 

id  of  suspended  sediment  size. 

4 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

5 

IDTYPT(NP,i) 

+ 

id  of  xy  series  to  describe  the  time -dependent  concentration 
profile  associated  with  the  Dirichlet  boundary  node  for  the 
suspended  sediment,  where  i  =  NCHEM+k. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

DB 

Card  group  identifier. 

C  3 

IC3 

3 

Dirichlet  boundary  conditions  for  particulate  chemicals  on 
suspended  sediments. 

1 

NP 

+ 

id  of  Dirichlet  boundary  node. 

2 

NPDBT(NP) 

+ 

id  of  global  node  corresponding  to  the  Dirichlet  boundary  node. 

3 

j 

+ 

id  of  suspended  sediment  size. 

4 

k 

+ 

id  of  particulate  chemical. 

5 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

6 

NSSIZE 

+ 

Number  of  sediment  sizes. 

7 

IDTYPT(NP,i) 

+ 

id  of  xy  series  to  describe  the  time -dependent  concentration  profile 
associated  with  the  Dirichlet  boundary  node  for  the  particulate 
chemical  on  the  suspended  sediment,  where  i  =  NCHEM  +  NSSIZE 
+  (j-1 )  *NCHEM  +  k,  j  6  [1, NSSIZE],  and  k  e[l,  NCHEM]. 

FB  Card 

Flux-tvoe 

boundary  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

FB 

Card  group  identifier. 

C  3 

IC3 

1 

Flux-type  boundary  conditions  for  dissolved  chemicals. 

1 

NS 

+ 

id  of  flux-type  boundary  side. 

2 

NSCBT(NS,1) 

+ 

id  of  global  node  corresponding  to  the  1-st  node  of  the  flux-type 
boundary  side. 

3 

NSCBT(NS,2) 

+ 

id  of  global  node  corresponding  to  the  2-nd  node  of  the  flux-type 
boundary  side. 

4 

NSCBT(NS,3) 

1 

Index  of  boundary  condition  type  for  the  flux-type  boundary  side. 
This  is  a  Variable  boundary  side  with  the  time-dependent 
concentration  specified  if  the  flow  is  incoming. 

2 

This  is  a  Cauchy  boundary  side  with  the  time-dependent  Cauchy 
flux  specified. 

3 

This  is  a  Neumann  boundary  side  with  the  time -dependent 
Neumann  flux  specified. 

5 

i 

+ 

id  of  dissolved  Chemical. 

6 

ICTYPT(NS,i) 

+ 

id  of  xy  series  to  describe  the  of  time-dependent  profile  associated 
with  the  flux-type  boundary  side  for  the  dissolved  chemical.  It  is 

(1)  a  concentration  profile  when  NSCBT(j,2)  =  1. 

(2)  a  Cauchy  flux  profile  when  NSCBT(j,2)  =  2. 

(3)  a  Neumann  flux  profile  when  NSCBT(j,2)  =  3. 

Field  Variable 

Value 

Description 

C  1-2  IC1 

FB 

Card  group  identifier. 
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C  3 

IC3 

2 

Flux-type  boundary  conditions  for  suspended  sediments. 

1 

NS 

+ 

id  of  flux-type  boundary  side. 

2 

NSCBT(NS,1) 

+ 

id  of  global  node  corresponding  to  the  1-st  node  of  the  flux-type 
boundary  side. 

3 

NSCBT(NS,2) 

+ 

id  of  global  node  corresponding  to  the  2-nd  node  of  the  flux-type 
boundary  side. 

4 

NSCBT(NS,3) 

1 

Index  of  boundary  condition  type  for  the  flux-type  boundary  side. 
This  is  a  Variable  boundary  node  with  the  time-dependent 
concentration  specified  if  the  flow  is  incoming. 

2 

This  is  a  Cauchy  boundary  node  with  the  time -dependent  Cauchy 
flux  specified. 

3 

This  is  a  Neumann  boundary  node  with  the  time-dependent 
Neumann  flux  specified. 

5 

k 

+ 

id  of  suspended  sediment  size. 

6 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

7 

ICTYPT(NS,i) 

+ 

id  of  xy  series  to  describe  the  time -dependent  profile  associated 

with  the  flux-type  boundary  side  for  the  suspended  sediment, 
where  i  =  NCHEM+k.  It  is 

(1)  a  concentration  profile  when  NSCBT(j,2)  =  1. 

(2)  a  Cauchy  flux  profile  when  NSCBT(j,2)  =  2. 

(3)  a  Neumann  flux  profile  when  NSCBT(j,2)  =  3. 


Field 

Variable 

Value 

Description 

C  1-2 

IC1 

FB 

Card  group  identifier. 

C  3 

IC3 

3 

Flux-type  boundary  conditions  for  particulate  chemicals  on 
suspended  sediments. 

1 

NS 

+ 

id  of  flux-type  boundary  side. 

2 

NSCBT(NS,1) 

+ 

id  of  global  node  corresponding  to  the  1-st  node  of  the  flux-type 
boundary  side. 

3 

NSCBT(NS,2) 

+ 

id  of  global  node  corresponding  to  the  2-nd  node  of  the  flux-type 
boundary  side. 

4 

NSCBT(NS,3) 

1 

Index  of  boundary  condition  type  for  the  flux-type  boundary  side. 
This  is  a  Variable  boundary  node  with  the  time -dependent 
concentration  specified  if  the  flow  is  incoming. 

2 

This  is  a  Cauchy  boundary  node  with  the  time -dependent  Cauchy 
flux  specified. 

3 

This  is  a  Neumann  boundary  node  with  the  time -dependent 
Neumann  flux  specified. 

5 

j 

+ 

id  of  suspended  sediment  size. 

6 

k 

+ 

id  of  particulate  chemical. 

7 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

8 

NSSIZE 

+ 

Number  of  sediment  sizes. 

9 

ICTYPT(NS,i) 

+ 

id  of  xy  series  to  describe  the  time -dependent  profile  associated 

with  the  flux-type  boundary  side  for  the  particulate  chemical  on  the 
suspended  sediment,  where  i  =  NCHEM  +  NSSIZE  +  (j- 
1)*NCHEM  +  k,  j  £ [1, NSSIZE],  and  k  £[1,  NCHEM].  It  is 

(1)  a  concentration  profile  when  NSCBT(j,2)  =  1. 

(2)  a  Cauchy  flux  profile  when  NSCBT(j,2)  =  2. 


1-14 


(3)  a  Neumann  flux  profile  when  NSCBT(j,2)  =  3. 


ON  Card 

Ocean- 

-type  boundary  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

ON 

Card  group  identifier. 

C  3 

IC3 

1 

Material  property  for  Ocean-type  boundary  conditions. 

1 

i 

+ 

id  of  material  type. 

2 

FDCAY(i) 

+ 

Flushing  decay  rate.  [1/T] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

ON 

Card  group  identifier. 

C  3 

IC3 

2 

Ocean-type  boundary  conditions  for  dissolved  chemicals. 

1 

NS 

+ 

id  of  flux-type  boundary  side. 

2 

NSNBT(NS,1) 

+ 

id  of  global  node  corresponding  to  the  1-st  node  of  the  ocean-type 
boundary  side. 

3 

NSNBT(NS,2) 

+ 

id  of  global  node  corresponding  to  the  2-nd  node  of  the  ocean-type 
boundary  side. 

4 

i 

+ 

id  of  dissolved  Chemical. 

5 

INTYPT(NS,i) 

+ 

material  id  associated  with  the  ocean-type  boundary  side  for 
specifying  the  flushing  decay  rate  of  the  dissolved  chemical. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

ON 

Card  group  identifier. 

C  3 

IC3 

3 

Ocean-type  boundary  conditions  for  suspended  sediments. 

1 

NS 

+ 

id  of  ocean-type  boundary  side. 

2 

NSNBT(NS,1) 

+ 

id  of  global  node  corresponding  to  the  1-st  node  of  the  ocean-type 
boundary  side. 

3 

NSNBT(NS,2) 

+ 

id  of  global  node  corresponding  to  the  2-nd  node  of  the  ocean-type 
boundary  side. 

4 

k 

+ 

id  of  suspended  sediment  size. 

5 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

6 

INTYPT(NS,i) 

+ 

material  id  of  associated  with  the  ocean-type  boundary  side  for 
specifying  the  flusing  decay  rate  of  the  suspended  sediment,  where 
i  =  NCHEM+k. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

ON 

Card  group  identifier. 

C  3 

IC3 

4 

Ocean-type  boundary  conditions  for  particulate  chemicals  on 
suspended  sediments. 

1 

NS 

+ 

id  of  ocean-type  boundary  side. 

2 

NSNBT(NS,1) 

+ 

id  of  global  node  corresponding  to  the  1-st  node  of  the  ocean-type 
boundary  side. 

3 

NSNBT(NS,2) 

+ 

id  of  global  node  corresponding  to  the  2-nd  node  of  the  ocean-type 
boundary  side. 

4 

j 

+ 

id  of  suspended  sediment  size. 

5 

k 

+ 

id  of  particulate  chemical. 

6 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

7 

NSSIZE 

+ 

Number  of  sediment  sizes. 
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8 

INTYPT(NS,i) 

+ 

material  id  of  associated  with  the  ocean-type  boundary  side  for 
specifying  the  flushing  decay  rate  of  the  particulate  chemical  on 
the  suspended  sediment,  where  i  =  NCHEM  +  NSSIZE  +  (j- 
1)*NCHEM  +  k,  j  6 [1, NSSIZE],  and  k  e[l,  NCHEM]. 

RS  card 

Rainfall  source  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

RS 

Card  group  identifier. 

C  3 

IC3 

1 

Rainfall  source  profile  for  dissolved  chemicals. 

1 

j 

+ 

id  of  dissolved  chemical. 

2 

lRC(j) 

+ 

id  of  xy  series  to  describe  the  time -dependent  rainfall  source 
profile  for  the  dissolved  chemical. 

IC  card 

Initial  condition  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IC 

Card  group  identifier. 

C  3 

IC3 

1 

Initial  condition  parameters. 

1 

IDCONC 

0 

index  of  initial  concentrations, 
constant  initial  concentrations. 

1 

variable  initial  concentrations.  In  this  case,  the  user  needs  to 
specify  ICCN,  ICCP,  and  ICCS  files  in  the  super  file  as  explained 
in  the  end  of  this  appendix. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IC 

Card  group  identifier. 

C  3 

IC3 

2 

Initially  constant  conditions  for  dissolved  chemicals. 

1 

j 

+ 

id  of  the  dissolved  chemical  being  input. 

2 

CCONST(j) 

+ 

Specified  initial  concentration  for  the  dissolved  chemical.  [M/L3] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IC 

Card  group  identifier. 

C3 

IC3 

3 

Initially  constant  conditions  for  particulate  chemicals  on  suspended 
sediments 

1 

j 

+ 

id  of  the  particulate  chemical  being  input. 

2 

i 

+ 

id  of  the  suspended  sediment  size  where  the  particulate  chemical 
is. 

Number  of  dissolved  chemicals. 

3 

NCHEM 

+ 

4 

CCONST(k) 

+ 

Specified  initial  concentration  for  the  particulate  chemical  on 
suspended  sediment  being  input,  where  k  =  i*NCHEM  +  j.  [M/M] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IC 

Card  group  identifier. 

C  3 

IC3 

4 

Initially  constant  conditions  for  particulate  chemicals  on  bed 
sediments 

1 

j 

+ 

id  of  the  particulate  chemical  being  input. 

2 

i 

+ 

id  of  the  bed  sediment  size  where  the  particulate  chemical  is. 

3 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

4 

NSSIZE 

+ 

Number  of  sediment  sizes. 
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5 

CCONST(k) 

+ 

Specified  initial  concentration  for  the  particulate  chemical  on 
the  bed  sediment  being  input,  where  k  =  (I+NSSIZE)* 

NCHEM +j.  [M/M] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IC 

Card  group  identifier. 

C  3 

IC3 

5 

Initially  constant  conditions  for  suspended  sediments. 

1 

i 

+ 

id  of  the  suspended  sediment  size  being  input. 

2 

SCONST(i) 

+ 

Specified  initial  concentration  for  the  suspended  sediment  of 
the  i-th  size.  [M/L3] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IC 

Card  group  identifier. 

C  3 

IC3 

6 

Initially  constant  conditions  for  bed  sediments. 

1 

i 

+ 

id  of  the  bed  sediment  size  being  input. 

2 

NSSIZE 

+ 

Number  of  sediment  sizes. 

3 

SCONST(k) 

+ 

Specified  initial  concentration  for  the  bed  sediment  of  the  i-th  size, 
where  k  =  NSSIZE  +  i.  [M/L3] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

IC 

Card  group  identifier. 

C  3 

IC3 

7 

Initially  constant  conditions  for  dissolved  chemicals  in  the 
interstitial  water  of  bed  sediments 

1 

i 

+ 

id  of  the  dissolved  chemical  being  input. 

2 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

3 

NSSIZE 

+ 

Number  of  sediment  sizes. 

4 

CCONST(k) 

+ 

Specified  initial  concentration  for  the  bed  sediment  of  the  i-th  size, 
where  k  =  NSSIZE  +  i.  [M/L3] 

XY  Card 

X-Y  Series  Parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

XY 

Card  group  identifier. 

C  3 

IC3 

1 

Generation  of  any  x-y  series  function. 

1 

i 

+ 

Index  number  for  x-y  series. 

2 

NPOINT(i) 

+ 

The  number  of  x-y  values  in  the  series. 

3 

DELTAX 

0 

Dummy  values. 

4 

DELTAY 

0 

Dummy  values. 

5 

REPEAT 

0 

Dummy  values. 

6 

BEGCYC 

0 

Dummy  values. 

7 

TNAME 

+ 

A  character  string  representing  the  name  of  the  XY  series. 

(new  line)X,Y 

After  XY  card,  the  x-y  values  of  the  series  are  listed  one  pair  per  line  up 

to  NPOINT(i). 

EN  Card 

End  of  data  control 

Field 

Variable 

Value 

Description 

Cl -2 

IC1 

EN 

Card  group  identifier. 

C  3 

IC3 

D 

End  of  input  data. 
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1.4.  File  4:  contaminant/sediment  and  reaction  information  file  (CHEM  file) 

T1-T3  Cards  Job  title 

Only  one  Tl,  T2,  and  T3  card  can  be  used. 

RX  card 
Field  Variable 
C  1-2  IC1 
C  3  IC3 

1  NRXN 

2  NCHEM 


Chemical  reaction  parameters 
Value  Description 

RX  Card  group  identifier. 

1  Numbers  of  reactions  and  chemicals 

+  Number  of  reactions. 

+  Number  of  dissolved  chemicals. 


Field  Variable  Value  Description 

C  1-2  IC1  RX  Card  group  identifier. 

C  3  IC3  2  Chemical  names 

1  NCHEM  +  Number  of  dissolved  Chemical. 

Note:  NCHEM  new  lines  are  needed. 

(the  first  new  line) 

CHEMN(l)  +  Name  of  the  1-st  dissolved  chemical, 

(the  second  new  line) 

CHEMN(2)  +  Name  of  the  2-nd  dissolved  chemical. 


(the  NCHEM-th  new  line) 

CHEMN(NCHEM)  +  Name  of  the  NCHEM-th  dissolved  chemical. 

Field  Variable  Value  Description 

C  1-2  IC1  RX  Card  group  identifier. 

C  3  IC3  3  Reaction  rate  constants  and  stoichiometry 

1  j  +  id  of  chemical  reaction. 

2  NCHEM  +  Number  of  dissolved  chemicals. 

3  FKRX(j)  +  Forward  reaction  rate  constant  of  the  reaction,  [reaction  dependent] 

4  BKRX(j)  +  Backward  reaction  rate  constant  of  the  reaction,  [reaction 

dependent] 

Note:  Two  new  lines  are  needed. 

(the  first  new  line) 

1  NURTS(j,l)  +  Stoichiometric  coefficient  of  the  1 -st  dissolved  chemical  (on  the 

reactant  side)  in  the  reaction. 

2  NURTS(j,2)  +  Stoichiometric  coefficient  of  the  2-nd  dissolved  chemical  (on  the 

reactant  side)  in  the  reaction. 


up  to  NURTS(j, NCHEM) 

(the  second  new  line) 

1  NUPDS(j,l)  +  Stoichiometric  coefficient  of  the  1 -st  dissolved  chemical  (on  the 

product  side)  in  the  reaction. 

2  NUPDS(j,2)  +  Stoichiometric  coefficient  of  the  2-nd  dissolved  chemical  (on  the 

product  side)  in  the  reaction. 
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up  to  NUPDS(j  ,NCHEM) 

ST  card  Sediment  parameters 


Field 

Variable 

Value 

Description 

C  1-2 

IC1 

ST 

Card  group  identifier. 

C  3 

IC3 

1 

Number  of  sediment  sizes 

1 

NSSIZE 

+ 

Number  of  sediment  sizes. 

2 

EXCHNG 

+ 

Exchange  coefficient  accounting  for  diffusion  of  chemicals 
between  pore  water  and  column  water. 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

ST 

Card  group  identifier. 

C  3 

IC3 

2 

Parameter  associated  with  each  sediment  size. 

1 

j 

+ 

id  of  Sediment  size. 

2 

SPARA(j,l) 

1 

The  type  of  the  sediment  size. 

Cohesive  sediment,  e.g.,  clay  and  silt. 

2 

Non-cohesive  sediment,  e.g.,  sand. 

3 

SPARA(j,2) 

+ 

Settling  speed  of  the  sediment  size.  [L/T] 

4 

SPARA(j,3) 

+ 

Critical  shear  stress  for  deposition  of  the  sediment  size.  [M/L/T2] 

5 

SPARA(j,4) 

+ 

Critical  shear  stress  for  erosion  of  the  j-th  size  fraction  of 
sediment.  [M/L/T2] 

6 

SPARA(j,5) 

+ 

Erodibility  of  the  sediment  size.  [M/L2] 

7 

SPARA(j,6) 

+ 

Specific  weight  of  the  sediment  size.  [M/L3] 

8 

SPARA(j,7) 

+ 

Median  diameter  of  particles  in  mixture  of  the  sediment  size.  [L] 

9 

SPARA(j,8) 

+ 

Critical  bottom  shear  stress  of  the  sediment  size  at  which  sediment 
movement  begins.  [M/L/T2] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

ST 

Card  group  identifier. 

C  3 

IC3 

3 

Sorption  rate  constants  related  to  suspended  sediment. 

1 

j 

+ 

id  of  suspended  sediment  size. 

2 

i 

+ 

id  of  dissolved  chemical. 

3 

FKSPSS(j,l) 

+ 

Forward  rate  constant  of  the  sorption  of  i-th  dissolved  chemical 
onto  the  suspended  sediment  of  the  j-th  size.  [L3/M/T] 

4 

BKSPSS(j,2) 

+ 

Backward  rate  constant  of  the  sorption  of  i-th  dissolved  chemical 
onto  the  suspended  sediment  of  the  j-th  size.  [1/T] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

ST 

Card  group  identifier. 

C  3 

IC3 

4 

Soiption  rate  constants  related  to  bed 
sediment. 

1 

j 

+ 

id  of  bed  sediment  size. 

2 

i 

+ 

id  of  dissolved  chemical. 

3 

FKSPBS(j,i) 

+ 

Forward  rate  constant  of  the  sorption  of  i-th  dissolved  chemical 
onto  the  bed  sediment  of  the  j-th  size.  [L3/M/T] 

4 

BKSPBS(j,i) 

+ 

Backward  rate  constant  of  the  sorption  of  i-th  dissolved  chemical 
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onto  the  bed  sediment  of  the  j-th  size.  [1/T] 


5. 

FKSP2(j,i) 

+ 

Forward  rate  constant  of  the  sorption  of  i-th  chemical  in  the 
interstitial  water  onto  the  bed  sediment  of  the  j-th  size.  [L3/M/T] 

6. 

BKSP2(j,i) 

+ 

Backward  rate  constant  of  the  sorption  of  i-th  chemical  in  the 
interstitial  water  onto  the  bed  sediment  of  the  j-th  size.  [L3/M/T] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

ST 

Card  group  identifier. 

C  3 

IC3 

5 

The  intrinsic  properties  related  to  bed  sediment. 

1 

j 

+ 

id  of  bed  sediment  size. 

2 

RHOB(j) 

+ 

Bulk  density  of  the  j-th  size  particle  of  bed  sediments.  [M/L3] 

3 

THETA(j) 

+ 

Porosity  of  the  j-th  size  particle  of  bed  sediments. 

VO  card 

Volatilization  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

VO 

Card  group  identifier. 

C  3 

IC3 

1 

Volatilization  rate  constants. 

1 

j 

+ 

id  of  dissolved  chemical. 

2 

IAP(j) 

+ 

id  of  xy  series  to  describe  the  time-dependent  partial  atmospheric 
pressure  profile  of  the  dissolved  chemical,  [atm] 

3 

FKVO(j) 

+ 

Forward  volatilization  rate  constant  of  the  dissolved  chemical. 
[1/T] 

4 

BKVO(j) 

+ 

Backward  volatilization  rate  constant  of  the  dissolved  chemical. 

[M/atm/L3/T] 


DY  card 

1  -st  order 

decay  parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

DY 

Card  group  identifier. 

C  3 

IC3 

1 

Decay  constants  of  dissolved  chemicals. 

1 

j 

+ 

id  of  dissolved  chemical. 

2 

DECAY(j) 

+ 

1-st  order  decay  constant  of  the  dissolved  chemical.  [1/T] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

DY 

Card  group  identifier. 

C  3 

IC3 

2 

Decay  constants  of  particulate  chemicals  on  suspended  sediments. 

1 

j 

+ 

id  of  sediment  size. 

2 

i 

+ 

id  of  particulate  chemical. 

3 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

4 

DCY 

+ 

Decay  constant  of  the  particulate  chemical  on  the  suspended 
sediment.  [1/T] 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

DY 

Card  group  identifier. 

C  3 

IC3 

3 

Decay  constants  of  particulate  chemicals  on  bed  sediments. 

1 

j 

+ 

id  of  sediment  size. 

2 

i 

+ 

id  of  particulate  chemical. 

3 

NCHEM 

+ 

Number  of  dissolved  chemicals. 

4 

NSSIZE 

+ 

Number  of  sediment  sizes. 
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5 

DCY 

+ 

Decay  constant  of  the  particulate  chemical  on  the  bed  sediment. 
[1/T] 

XY  Card 

X-Y  Series  Parameters 

Field 

Variable 

Value 

Description 

C  1-2 

IC1 

XY 

Card  group  identifier. 

C  3 

IC3 

1 

Generation  of  any  x-y  series  function. 

1 

i 

+ 

Index  number  for  x-y  series. 

2 

NPOINT(i) 

+ 

The  number  of  x-y  values  in  the  series. 

3 

DELTAX 

0 

Dummy  values. 

4 

DELTAY 

0 

Dummy  values. 

5 

REPEAT 

0 

Dummy  values. 

6 

BEGCYC 

0 

Dummy  values. 

7 

TNAME 

+ 

A  character  string  representing  the  name  of  the  XY  series. 

(new  line)X,Y 

After  XY  card,  the  x-y  values  of  the  series  are  listed  one  pair  per  line  up 

to  NPOINT(i). 

EN  Card 

End  of  data  control 

Field 

Variable 

Value 

Description 

Cl -2 

IC1 

EN 

Card  group  identifier. 

C  3 

IC3 

D 

End  of  input  data. 

1.5.  File  5:  initial  water  surface  elevation  file  (ICWD  file) 

This  file  includes  the  following  lines: 

1st  line:  (A7) 

1 .  File  title. 

2nd  line:  (A7,  IX,  A6) 

1 .  Data  title. 

2.  BEGSCL 

3rd  line:  (A6) 

1  BEGSCL 

4th  line:  (A2.I8) 

1  ND 

2.  number  of  nodes. 

5th  line:  (A2,I8) 

1.  EM 

2.  number  of  elements. 
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6th  line:  (A4,1X,A40) 

1.  subtitle  1. 

2.  subtitle  2. 

7th  line:  (A2,A4,D16.8) 

1  TI 

2.  any  four  characters. 

3.  initial  time. 

In  the  next  NNP  (number  of  nodes)  lines,  each  line  contains  one  record  stating  the  water  surface  elevation 
at  the  corresponding  node. 

1 .  initial  water  surface  elevation. 


1.6.  File  6:  initial  velocity  (ICQL  file) 

This  file  includes  the  following  lines: 

1st  line:  (A7) 

2.  File  title. 

2nd  line:  (A7,  IX,  A6) 

1 .  Data  title. 

2.  BEGVEC 

3rd  line:  (A6) 

1  BEGVEC 

4th  line:  (A2.I8) 

1  ND 

2.  number  of  nodes. 

5th  line:  (A2,I8) 

1.  EM 

2.  number  of  elements. 

6th  line:  (A4,1X,A40) 

1.  subtitle  1. 

2.  subtitle  2. 

7th  line:  (A2,A4,D16.8) 

1  TI 

2.  any  four  characters. 

3.  initial  time. 


In  the  next  NNP  (number  of  nodes)  lines,  each  line  contains  two  record  stating  the  x-  and  y-  velocity 
components  at  the  corresponding  node.  (Free  format) 


1. 

2. 
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initial  x-velocity. 
initial  y-velocity. 


1.7.  File  7:  initial  dissolved  chemical  concentration  file  (ICCN  file) 

This  file  includes  the  following  lines 

1st  line:  (A7) 

3.  File  title. 

2nd  line:  (A7,  IX,  A6) 

1 .  Data  title. 

2.  BEGSCL 

3rd  line:  (A6) 

1  BEGSCL 

4th  line:  (A2,I8) 

1  ND 

2.  number  of  nodes. 

5th  line:  (A2.I8) 

1.  EM 

2.  number  of  elements. 

6th  line:  (A4,1X,A40) 

1.  subtitle  1. 

2.  subtitle  2. 

7th  line:  (A2,A4,D16.8) 

1  TI 

2.  any  four  characters. 

3.  initial  time. 

In  the  next  NCHEM*NNP  (number  of  nodes)  lines,  the  user  input  the  dissolved  concentrations  according 
to  the  following  logic.  (Free  format) 

DO  I  =  1,NCHEM 
DO  NP=1,NNP 
C(NP,I)  CB(NP,I) 

ENDDO 

ENDDO 

where  C(NP,  I)  and  CB(NP,I)  are  the  initial  dissolved  concentration  in  the  water  column  and  in  the  interstitial 
water,  respectively,  associated  with  the  Ilh  chemical  at  the  NPth  node. 


1.8.  File  8:  initial  particulate  chemical  concentration  file  (ICPN  file) 
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This  file  includes  the  following  lines 

1st  line:  (A7) 

1 .  File  title. 

2nd  line:  (A7,  IX,  A6) 

1 .  Data  title. 

2.  BEGSCL 

3rd  line:  (A6) 

1  BEGSCL 

4th  line:  (A2,I8) 

1  ND 

2.  number  of  nodes. 

5th  line:  (A2.I8) 

1.  EM 

2.  number  of  elements. 

6th  line:  (A4,1X,A40) 

1.  subtitle  1. 

2.  subtitle  2. 

7th  line:  (A2,A4,D16.8) 

1  TI 

2.  any  four  characters. 

3.  initial  time. 

In  the  next  NCHEM*NSSIZE*NNP  (number  of  nodes)  lines,  the  user  input  the  particulate  chemical 
concentrations  according  to  the  following  logic.  (Free  format) 

DO  I  =  1,NCHEM 
DO  J=1,NSSIZE 
DO  NP=1,NNP 
PS(NP,J,I)  PB(NP,J,I) 

ENDDO 

ENDDO 

ENDDO 

where  PS(NP,J,I)  and  PB(NP,J,I)  are  the  initial  particulate  chemical  concentration  associated  with  the  Ith 
chemical  on  the  J'h  suspended  and  bed  sediments,  respectively,  at  the  NPth  node. 


1.9.  File  9:  initial  sediment  concentration  file  (ICSN  file) 

This  file  includes  the  following  lines 
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1st  line:  (A7) 

1 .  File  title. 

2nd  line:  (A7,  IX,  A6) 

1 .  Data  title. 

2.  BEGSCL 

3rd  line:  (A6) 

1  BEGSCL 

4th  line:  (A2,I8) 

1.  ND 

2.  number  of  nodes. 

5th  line:  (A2.I8) 

1.  EM 

2.  number  of  elements. 

6th  line:  (A4,1X,A40) 

1.  subtitle  1. 

2.  subtitle  2. 

7th  line:  (A2,A4,D16.8) 

1  TI 

2.  any  four  characters. 

3.  initial  time. 

In  the  next  NSSIZE*NNP  (number  of  nodes)  lines,  the  user  input  the  sediment  concentrations  according  to 
the  following  logic.  (Free  format) 

DO  J=1,NSSIZE 
DO  NP=1,NNP 
SS(NP,J)  BS(NP,J) 

ENDDO 

ENDDO 

where  SS(NP,J)  and  BS(NP,J)  are  the  initial  suspended  and  bed  sediment  concentrations,  respectively,  at  the 
NPlh  node. 

In  the  table  below,  we  give  brief  description  on  the  files  that  should  be  included  in  the  super  files 
for  achieving  desired  simulations.  In  the  super  file,  the  BCFT  file  is  a  MUST  in  order  to  read  in  the 
simulation  control  index,  KMOD,  as  defined  above.  Both  the  BCTT  and  CHEM  files  need  to  be  included 
in  the  super  file  when  transport  simulations  are  desired.  When  variable  initial  conditions  are  to  be  read,  more 
files  need  to  be  taken  into  account  in  the  super  files,  as  stated  in  the  table  below. 
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Index  in  the 
super  file 

Data  File  Description 

GEOM 

geometry  file 

BCFT 

flow  file 

BCTT 

transport  file 

CHEM 

Chemistry  file 

ICWD 

initial  water  surface  elevation 

ICQL 

initial  velocity 

ICCN 

initial  concentrations  of  dissolved  chemicals:  data  are  read  chemical  by  chemical;  for 
each  chemical,  the  concentration  in  the  water  column  is  listed  first,  followed  by  that  in 
the  interstitial  water 

ICPN 

initial  concentrations  of  particulate  chemicals:  data  are  read  chemical  by  chemical;  for 
each  particulate  chemical  on  a  sediment  size,  the  particulate  chemical  concentration  on 
the  suspended  sediment  is  first  listed,  then  the  corresponding  one  on  the  bed  sediment 

ICSN 

initial  concentrations  of  sediments:  data  are  read  size  by  size;  for  each  sediment  size,  the 
suspended  sediment  concentration  is  first  listed,  followed  the  bed  sediment  concentration 

HSOL 

computational  results  of  flow,  including  both  water  surface  elevation,  velocity,  and  depth 

CSOL 

computational  results  of  dissolved  chemical  concentrations:  for  each  chemical,  the 
concentration  in  the  water  column  is  listed  first,  followed  by  that  in  the  interstitial  water 

PSOL 

computational  results  of  particulate  chemical  concentrations:  for  each  particulate 
chemical  on  a  sediment  size,  the  particulate  chemical  concentration  on  the  suspended 
sediment  is  first  listed,  then  the  corresponding  one  on  the  bed  sediment 

SSOL 

computational  results  of  sediment  concentrations:  for  each  sediment  size,  the  suspended 
sediment  concentration  is  first  listed,  followed  the  bed  sediment  concentration 

Note:  All  files  above  are  ASCII  files.  We  also  store  computational  results  in  a  binary  file,  best.sto,  for 
plotting  purposes. 


Here  we  provide  an  example  to  show  the  contents  of  a  super  file.  Suppose  we  are  to  simulate  both 

flow  and  transport,  the  initial  conditions  are  variable,  and  we  want  to  print  the  computational  results  of  flow 

as  well  as  transport,  then  we’ll  have  a  super  file  as  follows. 

WMS2DSUP 
GEOM  test.2dm 
BCFT  test.2bc 
BCTT  test.2tp 
CHEM  test.2ch 
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ICWD 

test.ich 

ICQL 

test.icv 

ICCN 

test.icc 

ICPN 

test.icp 

ICSN 

test.ics 

HSOL 

test_f.dat 

CSOL 

test_c.dat 

PSOL 

test_p.dat 

SSOL 

test_s.dat 

Thus,  we  have  test.2dm  describing  the  geometry,  test.2bc  specify  the  flow  simulation,  test.2tp  specify 
the  transport  simulation,  test.2ch  provide  contaminant/sediment  as  well  as  reaction  information,  test.ich, 
test.icv,  test.icc,  test.icp,  and  test.ics  provide  the  initial  conditions,  and  test_f.dat,  test_c.dat,  test_p.dat,  and 
test_s.dat  record  the  computational  results  at  desire  times.  The  following  give  the  idea  of  how  computational 
results  are  recorded  in  HSOL,  CSOL,  PSOL,  and  SSOL  files. 


HSOL  file:  hydrodynamic  variables 

DO  NP  =  1,NNP 

NP  X(NP)  Y(NP)  ETA(NP)  U(NP)  V(NP)  D(NP) 
ENDDO 


CSOL  file:  Dissolved  chemicals  in  the  water  column  and  the  interstitial  water 

DO  I  =  1,NCHEM 
DO  NP=1,NNP 

I  NP  X(NP)  Y(NP)  C(NP,I)  CB(NP,I) 

ENDDO 

ENDDO 


PSOL  file:  Particulate  chemicals  on  suspended  and  bed  sediments 

DO  I  =  1,NCHEM 
DO  J=1,NSSIZE 
DO  NP=1,NNP 

I  J  NP  X(NP)  Y(NP)  PS(NP,J,I)  PB(NP,J,I) 

ENDDO 

ENDDO 

ENDDO 


SSOL  file:  Suspended  and  bed  sediments 

DO  J=1,NSSIZE 
DO  NP=1,NNP 

J  NP  X(NP)  Y(NP)  SS(NP,J)  BS(NP,J) 
ENDDO 
ENDDO 
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Appendix  II.  PROGRAM  STRUCTURE  AND  SUBROUTINE  DESCRIPTION 

In  this  appendix,  we  are  to  give  plots  showing  the  structure  of  the  computer  program  and 
provide  description  for  each  subroutine. 

II.l.  Program  structure 

1.  Main  program  and  subroutines  it  calls  to. 


II-2 


2.  Subroutines  concerning  the  computation  of  water  flow. 


3.  Subroutines  concerning  the  computation  of  chemical/sediment  transport. 
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4.  Subroutines  concerning  backward  particle  tracking. 


5.  Subroutines  concerning  particle  tracking  in  a  subelement. 
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II.2.  Subroutine  description 
MAIN 

This  main  program  serves  the  purpose  of  controlling  the  flow  path  of  desired  simulations.  It  calls 
to  Subroutine  FILES_2  to  read  in  necessary  data  to  undergo  numerical  simulations.  It  calls  to  Subroutine 
FLOW2MOC  to  compute  2-D  flow  and  Subroutine  TRANSP_2  to  calculate  2-D  sediment  and  chemical 
transprot.  In  addition,  it  calls  to  needed  subroutines  to  print  or  store  numerical  results  at  desired  times. 

Subroutine  ACTEST 

This  subroutine  is  used  to  determine  a  quasi-steady  flow  by  checking  the  maximum  absolute  error 
of  water  stages  between  two  successive  tidal  cycles. 

Subroutine  ADVBCF2M 

This  subroutine  is  to  implement  boundary  conditions  for  the  computation  in  the  Lagrangian  step  of 
solving  the  2-D  flow  equations  when  the  method  of  characteristics  is  applied.  In  this  subroutine,  boundary 
conditions  with  either  water  stage  or  inflow  flux  specified  are  implemented  to  obtain  water  stage  or  normal 
flux,  respectively,  at  associated  boundary  nodes. 

Subroutine  ADVBCT_2 

This  subroutine  is  to  implement  boundary  conditions  for  the  computation  in  the  Lagrangian  step  of 
solving  the  2-D  transport  equations  of  (1)  suspended  sediment,  (2)  dissolved  chemicals,  and  (3)  particulate 
chemicals  on  suspended  sediments.  Since  the  diffusive  flux  is  neglected  in  the  Lagrangian  Step,  only  are 
Dirichlet,  Cauchy,  and  variable  boundary  conditions  are  applied. 

Subroutine  ADWK_2 

This  subroutine  is  to  determine  the  working  arrays  of  coordinates  (XW),  velocities  (VXW), 
subelement  indices  (IEW),  subelement  boundary  index  (1BW),  and  subelement  connectivity  (NLRLW  and 
LRLW)  in  the  element  being  considered  in  order  to  perform  the  2-D  “in-element”  particle  tracking. 

Subroutine  ALGBY_2 

This  subroutine  is  to  achieve  particle  tracking  along  the  2-D  unspecified  boundary.  Water  flow  is 
considered  parallel  to  the  boundary  on  the  unspecified  boundary.  This  subroutine  also  considers  the  first- 
order  growth/decay  term  along  characteristics. 

Subroutine  ALGBYF2M 

This  subroutine  is  to  achieve  particle  tracking  along  the  2-D  unspecified  boundary.  Water  flow  is 
considered  parallel  to  the  boundary  on  the  unspecified  boundary.  It  is  used  only  when  the  method  of 
characteristics  is  employed  to  solve  the  2-D  flow  equation. 

Subroutine  AREA 

This  subroutine  is  to  compute  the  area  of  a  2-D  element,  either  quadrilateral  or  triangular. 

Subroutine  ASEMBF2M 

This  subroutine  is  to  compose  the  global  matrix  equation  to  adjust  2-D  velocity  in  the  last  part  of  the 
method  of  characteristic  approach  when  eddy  viscosity  is  taken  into  account. 
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Subroutine  ASEMBT_2 

This  subroutine  is  to  compute  global  coefficient  matrices  CMATRX1  and  CMATRX2  for  2-D 
transport  equations.  CMATRX1  is  a  mass  matrices  that  takes  into  account  the  total  time  derivative  term, 
while  CMATRX2  is  stiffness  matrices  that  considers  the  dispersion  term. 

Subroutine  ASOLV2M 

This  subroutine  is  to  construct  and  solve  a  set  of  three  algebraic  equations  (Eqs.  (3.129)  through 
(3.131))  that  is  formed  by  using  the  method  of  characteristics  to  solve  the  dynamic  wave  model. 

Subroutine  BASE_21 

This  subroutine  is  to  compute  the  shape  functions  and  their  derivatives  at  a  point  located  in  a 
quadrilateral  element,  given  the  local  coordinates  at  the  point. 

Subroutine  BASE_22 

This  subroutine  is  to  compute  the  shape  functions  at  a  point  located  in  a  2-D  element  (either 
triangular  or  quadrilateral),  given  the  Cartesian  coordinates  both  at  the  point  and  element  nodes. 

Subroutine  BCT_2 

This  subroutine  is  to  implement  boundary  conditions,  including  Dirichlet,  Cauchy,  Neumann,  and 
variable  types,  for  2-D  transport.  Dirichlet  boundary  conditions  are  specified  at  boundary  nodes,  while  the 
other  types  of  boundary  conditions  are  categorized  as  flux-type  boundary  conditions  and  are  specified  on 
boundary  sides.  Dirichlet,  Cauchy,  Neumann,  and  conventional  variable  boundary  conditions  are  prescribed 
with  time -dependent  profiles. 

Subroutine  BSHEAR_2 

This  subroutine  is  to  achieve  the  calculation  of  bottom  shear  stress  at  all  nodes. 

Subroutine  BTRAKF2M 

This  subroutine  is  to  implement  backward  particle  tracking  to  obtain  the  water  depth  at  all  2-D  global 
nodes.  Since  a  first-order  approach  is  used  for  both  kinematic  and  diffusion  models  and  all  source/sink  rate 
terms  can  be  expressed  a  linear  function  of  water  depth,  we  thus  take  into  account  the  first-order 
decay/growth  rate  along  the  characteristics  when  we  implement  backward  particle  tracking.  As  for  the 
zeroth-order  source/sink  rate,  we  simply  multiply  this  rate  by  the  time-step  size  being  considered  and  add 
the  product  to  the  outcome  of  backward  particle  tracking. 

Subroutine  BTRAKT_2 

This  subroutine  is  to  implement  backward  particle  tracking  to  obtain  the  Lagrangian  values  of 
concentrations  of  (1)  dissolved  chemicals,  (2)  suspended  sediments,  and  (3)  particulate  chemicals  on 
suspended  sediments  at  all  2-D  global  nodes. 

Subroutine  CHNGSN 

This  subroutine  is  to  change  sign  of  the  velocities  at  point  P  (the  starting  location),  point  Q  (the 
ending  location),  and  working  subelement  nodes  in  the  process  of  2-D  or  3-D  “in-element”  particle  tracking. 

Subroutine  CKCNEL_2 

This  subroutine  is  to  determine  the  element  that  connects  to  a  specified  element  through  the  element 
side  with  the  Nl-th  and  N2-th  global  nodes  serving  as  the  two  end  nodes  of  the  2-D  side. 
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Subroutine  CKCOIN_2 

This  subroutine  is  to  determine  whether  or  not  a  specific  node  coincides  with  the  N 1  -th  global  node 
or  the  N2-th  global  node. 

Subroutine  CRACKD 

This  subroutine  is  to  read  real  numbers  in  order  from  a  line  data  record. 

Subroutine  CRACKI 

This  subroutine  is  to  read  integer  numbers  in  order  from  a  line  data  record. 

Subroutine  DEPERO_2 

This  subroutine  is  to  compute  the  deposition  and  erosion  rates  either  at  a  global  node.  Both  cohesive 
and  non-cohesive  cases  are  taken  into  account.  Eqs.  (3.129)  and  (3.130)  are  employed  for  cohesive 
sediments,  while  Eqs.  (3.131)  through  (3.134)  are  used  for  non-cohesive  sediments. 

Subroutine  DGELG 

This  subroutine  is  to  solve  matrix  equations  by  a  direction  solver  incorporated  with  full-pivoting. 

Subroutine  DISPC_2 

This  subroutine  is  to  compute  dispersion  coefficients  for  2-D  transport  simulations.  The  coefficient 
is  evaluated  at  every  quadrature  node  of  elements. 

Function  DOTPRD 

This  function  is  to  compute  the  inner  product  of  two  given  vectors. 

Subroutine  ELENOD_2 

This  subroutine  is  to  determine  the  number  of  element  nodes  for  a  specified  2-D  element. 

Subroutine  ELTRK_2 

This  subroutine  is  to  implement  the  2-D  “in-element”  particle  tracking  in  a  specified  element.  The 
first-order  growth/decay  term  is  also  taken  into  account  in  the  tracking. 

Subroutine  ELTRK2M 

This  subroutine  is  to  implement  the  2-D  “in-element”  particle  tracking  in  a  specified  element.  It  is 
used  only  when  the  method  of  characteristics  is  used  to  solve  the  2-D  flow  equation. 

Subroutine  EQGENT_2 

This  subroutine  is  to  generate  matrix  equations  for  2-D  transport  simulations.  When  IDC  =  0,  it  is 
for  the  transport  of  the  IDS-th  suspended  sediment.  When  IDS  =  0,  it  is  for  the  transport  of  the  IDC-th 
dissolved  chemical.  When  neither  IDC  nor  IDs  equals  zero,  it  is  for  the  IDC-th  particulate  chemical  on  the 
IDS-th  suspended  sediment. 

Subroutine  FBTABL 

This  subroutine  is  to  write  out  mass  balance  table  of  flow  at  desired  times. 

Function  FCOS_2 

This  is  to  determine  the  direction  of  the  outer  product  of  two  2-D  vectors. 
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Subroutine  FFLUX 

This  subroutine  is  to  compute  the  normal  flow  flux  through  each  boundary  node. 

Subroutine  FILES_2 

This  subroutine  is  to  read  in  information  in  order  to  perform  2-D  simulations.  The  information 
includes  (1)  geometry  data,  (2)  flow  simulation  data,  (3)  transport  simulation  data,  and  (4)  1-D/2-D  and  2- 
D/3-D  interface  data. 


Subroutine  FIXCKF2M 

This  subroutine  is  to  handle  the  case  in  the  particle  tracking  of  2-D  water  flow  when  a  particle 
encounters  the  boundary.  If  this  boundary  is  a  specified  (flow-through)  boundary,  the  current  particle 
tracking  stops  here  and  interpolation  is  employed  to  compute  the  Lagrangian  value.  If  this  boundary  is 
unspecified  (closed),  on  the  other  hand,  the  flow  direction  is  parallel  to  the  boundary  and  particle  tracking 
along  the  boundary  proceeds  until  either  the  available  tracking  time  is  completely  consumed  or  a  specified 
boundary  is  reached.  It  is  used  only  when  the  method  of  characteristics  is  employed  to  solve  the  2-D  flow 
equation. 

Subroutine  FIXCKT_2 

This  subroutine  is  to  handle  the  case  in  the  particle  tracking  of  2-D  chemical  transport  when  a 
particle  encounters  the  boundary.  If  this  boundary  is  a  specified  (flow-through)  boundary,  the  current 
particle  tracking  stops  here  and  interpolation  is  employed  to  compute  the  Lagrangian  value.  If  this  boundary 
is  unspecified  (closed),  on  the  other  hand,  the  flow  direction  is  parallel  to  the  boundary  and  particle  tracking 
along  the  boundary  proceeds  until  either  the  available  tracking  time  is  completely  consumed  or  a  specified 
boundary  is  reached. 


Function  FKAPA_2 

This  function  is  to  compute  the  bottom  friction  coefficient,  K,  as  defined  as  follows. 


K 


gn|v(h  +  Zo)|1/2|l+(VZo) 


2/3 


2/3 
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Subroutine  FLOW2MOC 

This  subroutine  serves  as  the  control  panel  of  computing  2-D  flow  by  using  the  method  of 
characteristics  to  solve  the  dynamic  model.  The  Picard  method  is  used  to  deal  with  nonlinearity.  During 
each  nonlinear  iteration,  The  characteristic  equation  (Eq.  (3 . 1 5))  is  solved  first  by  the  backward  particle 
tracking  along  three  characteristics.  After  H#,  u#,  and  v#  are  determined,  we  compute  H,  u,  and  v  by  solving 
Eqs.  (3.64)  through  (3.66).  Convergence  is  determined  by  checking  the  maximum  relative  error  of  H,  u,  and 
v.  Relaxation  technique  is  applied  to  update  the  working  water  stage  and  velocity  for  the  next  iteration  if 
convergence  has  not  been  reached. 

Subroutine  FLUX 

This  subroutine  is  to  compute  mass  flux  of  the  given  mobile  material  at  boundary  nodes. 

Subroutine  FSFLOW2 

This  subroutine  is  to  compute  the  data  needed  in  the  mass  balance  table  of  flow. 


Subroutine  FTAUB_2 
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This  subroutine  is  to  calculate  bottom  shear  stress  at  an  node.  The  bottom  shear  stress  is  computed 
according  to  the  following  equations. 


|V(h  +  Z„)| 

1/2 

1+(VZ0)2 

2/3  h  1/3 

(II- 2) 

< 
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N 

O 

1/2 

l+(VZo)2 

2/3  h  1/3 

(II- 3) 

where  Txb  and  Tyb  are  bottom  shear  stress  in  x-  and  y-directions,  respectively;  u  and  v  are  velocity  components 
in  x-  and  y-directions,  respectively;  p  is  water  density;  g  is  gravity;  n  is  Manning’s  n;  h  is  water  depth;  Z0 
is  bottom  elevation. 


Subroutine  GAUPNT_2 

This  subroutine  is  to  compute  shape  functions  and  their  derivatives  at  quadrature  nodes. 

Subroutine  GEOM_2 

This  subroutine  is  to  read  in  2-D  geometry,  including  coordinates  and  element  indices. 

Subroutine  GRAD_2 

This  subroutine  is  to  compute  the  spatial  derivatives  of  water  depth,  stage,  flow  velocity,  and 
Manning’s  n  as  desired  in  2-D  simulations.  In  addition,  the  second  spatial  derivatives  of  water  depth  and 
stage  are  calculated  as  needed. 

Subroutine  IBWCHK_2 

This  subroutine  is  to  determine  the  two  global  nodes  that  serve  as  the  two  end  node  of  an  element 
side.  This  element  side  is  the  side  that  the  particle  has  encountered  when  it  passes  through  the  element  being 
considered. 


Subroutine  IDFBS 

This  subroutine  is  to  relate  flux-type  boundary  sides  with  global  boundary  sides.  The  relationship 
is  stored  in  two  arrays:  one  for  flow  simulations  and  the  other  for  transport  simulations. 

Subroutine  ILUCG 

This  subroutine  is  to  solve  the  linearized  matrix  equation  that  is  sparse  asymmetric  with  the 
preconditioned  conjugate  gradient  method  using  the  incomplete  Cholesky  decomposition  as  a  preconditioner. 

Subroutine  INTRPF2M 

This  subroutine  is  to  determine  the  Lagrangian  characteristic  quantities  by  interpolation  after 
backward  particle  tracking  in  solving  flow  governing  equations  when  the  method  of  characteristics  is  used 
to  deal  with  the  dynamic  wave  model. 

Subroutine  INTRPT_2 

This  subroutine  is  to  determine  the  Lagrangian  chemical  and  suspended  sediment  concentrations  by 
interpolation  after  backward  particle  tracking  in  the  Lagrangian  step  of  solving  2-D  transport  governing 
equations. 


Subroutine  JACOBI_2 

This  subroutine  is  to  compute  the  Jacobian  for  the  chemical  system  at  the  NP-th  node.  The  Jacobian 


II-9 


elements  are  defined  as  in  Eqs.  (3.1 13)  through  (3.128). 

Subroutine  KXKY_2M 

This  subroutine  is  to  determine  characteristic  directions  at  a  2-D  global  node  for  flow  simulations, 
based  on  the  velocity,  the  spatial  velocity  derivatives,  and  the  spatial  stage  derivatives  at  the  node.  This 
subroutine  is  needed  only  when  the  method  of  characteristics  is  employed  to  solve  the  2-D  flow  equation. 

Subroutine  KXKYNP2M 

This  subroutine  is  to  determine  characteristic  directions  at  all  2-D  global  nodes  for  flow  simulations. 
This  subroutine  is  needed  only  when  the  method  of  characteristics  is  employed  to  solve  the  2-D  flow 
equation. 

Subroutine  LLTINV 

This  subroutine  is  to  solve  for  a  modified  residual  that  is  to  be  used  in  the  preconditioned  conjugate 
gradient  algorithm. 

Subroutine  LNDGEN_2 

This  subroutine  is  to  determine  node -node  and  node-element  connectivity  in  2-D  simulations,  based 
on  the  given  element  indices. 

Subroutine  LOCQ_2 

This  subroutine  is  to  locate  the  target  point  and  its  velocity  for  a  tracking  path  in  a  2-D  subelement 
of  particle  tracking.  The  available  tracking  time  left  after  this  tracking  path  is  also  computed. 

Subroutine  LOCQ2M 

This  subroutine  is  to  locate  the  target  point  and  its  velocity  for  a  tracking  path  in  a  2-D  subelement 
of  particle  tracking.  The  available  tracking  time  left  after  this  tracking  path  is  also  computed.  It  is  used  since 
the  method  of  characteristics  is  employed  to  solve  the  2-D  flow  equation. 

Subroutine  MMLOC_2 

This  subroutine  is  to  specify  the  starting  location  to  start  the  2-D  “in-element”  particle  tracking  in 
an  element. 

Subroutine  NELIM2 

This  subroutine  is  to  prepare  an  array  consisting  of  nodal  numbers  from  a  set  of  element  face  data. 

Subroutine  NODV_2 

This  subroutine  is  to  calculate  Manning’s  n  at  every  2-D  global  node. 

Subroutine  OBCOEF2M 

This  subroutine  is  to  determine  the  coefficient  of  a  characteristic  at  an  open  boundary  node.  The 
characteristic  has  been  determined  coming  from  the  outside  of  the  domain.  This  subroutine  is  needed  only 
when  the  method  of  characteristics  approach  is  taken  to  solve  2-D  flow  equations. 


Subroutine  ONLINE_2 

This  subroutine  is  to  adjust  the  coordinates  of  a  specified  point  if  needed,  such  that  this  point  will 
be  located  on  the  line  segment  with  two  given  2-D  nodes  serving  as  its  end  nodes. 
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Subroutine  PAGEN_2 

This  subroutine  is  to  generate  2-D  pointer  arrays  needed  for  numerical  simulations. 

Subroutine  PFFCT2 

This  subroutine  is  to  determine  the  functional  values  of  an  array  by  interpolation  from  a  series  of 
given  profiles. 

Subroutine  PISS 

This  subroutine  is  to  solver  the  linearized  matrix  equation  with  a  pointwise  iteration  solution 
strategy. 

Subroutine  POLYP 

This  subroutine  is  to  solve  a  modified  residual  that  is  to  be  used  in  the  preconditioned  conjugate 
gradient  algorithm. 

Subroutine  PPCG 

This  subroutine  is  to  solve  the  linearized  matrix  equation  with  the  preconditioned  conjugate  gradient 
method  by  using  a  polynomial  as  a  preconditioner. 

Subroutine  PREPAR2M 

This  subroutine  is  to  prepare  for  data  needed  to  construct  and  solve  the  set  of  algebraic  equations 
in  Subroutine  ASOLV2M. 

Subroutine  Q2FLUX 

This  subroutine  is  to  perform  boundary  integration  over  a  2-D  boundary  side. 

Subroutine  Q3D 

This  subroutine  is  to  compute  element  flux  matrices  due  to  dispersion  for  mobile  materials,  where 
the  element  is  triangular. 

Subroutine  Q3T 

This  subroutine  is  to  determine  the  triangular  element  coefficient  matrices  and  load  vectors  to  sole 
2-D  transport  equation  in  the  Eulerian  step. 

Subroutine  Q4D 

This  subroutine  is  to  compute  element  flux  matrices  due  to  dispersion  for  mobile  materials,  where 
the  element  is  quadrilateral. 

Subroutine  Q4T 

This  subroutine  is  to  determine  the  quadrilateral  element  coefficient  matrices  and  load  vectors  to  sole 
2-D  transport  equation  in  the  Eulerian  step. 
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Subroutine  Q34FM 

This  subroutine  is  to  determine  the  element  coefficient  matrices,  including  mass  matrices  and  the 
stiffness  matrices  taking  into  account  eddy  viscosity  in  soling  2-D  dynamic  wave  model. 

Subroutine  Q34F1 

This  subroutine  is  to  determine  element  coefficient  matrices  and  load  vectors  to  compute  the 
gradients  of  bottom  elevation,  water  stage,  velocities,  and  Manning’s  n  in  2-D  flow  simulations. 

Subroutine  Q34F2 

This  subroutine  is  to  determine  element  coefficient  matrices  and  load  vectors  to  compute  the 
gradients  of  the  derivatives  of  bottom  elevation  in  2-D  flow  simulations. 

Subroutine  RAINOD_2 

This  subroutine  is  to  determine  rainfall  rate  at  2-D  global  nodes. 

Subroutine  REPLAS_2 

This  subroutine  is  to  write  the  values  of  the  first  four  arguments  into  the  last  four  arguments. 

Subroutine  RESIDU_2 

This  subroutine  is  to  compute  the  residual  of  the  2-D  transport  equations  of  (1)  dissolved  chemicals 
(when  ID  =  1),  (2)  particulate  chemicals  on  suspended  sediments  (when  ID  =  1),  (3)  particulate  chemicals 
on  bed  sediments  (when  ID  =  1),  and  (4)  suspended  sediments  (when  ID  =  2)  at  the  NP-th  node.  The 
residuals  are  defined  as  in  Eqs.  (3.109)  through  (3.112)  for  contaminant  transport  equations  and  in  the 
following  equation  for  the  sediment  transport  equations. 

RESSn  =  h(Sn)N+1/2  +  At  (RHSSn)N+1  -  At  (RHSSn)N  ne  [1,NS]  (II.4) 


where  (RHSSn)N+1  is  evaluated  with  ID  =  22  in  Subroutine  RHSCOMP 1 ,  whereas  (RHSSn)N  is  calculated  with 
ID  =  21. 


Subroutine  RHSCMP  2 

This  subroutine  is  to  compute  the  right-hand  side  of  the  2-D  transport  governing  equations  of  (1) 
dissolved  chemicals  (when  ID  =  1),  (2)  particulate  chemicals  on  suspended  sediments  (when  ID  =  1),  (3) 
particulate  chemicals  on  bed  sediments  (when  ID  =  1),  and  (4)  suspended  sediments  (when  ID  =  21  or  22) 
at  the  NP-th  node.  The  right-hand  sides  are  defined  as  in  Eqs.  (3.102),  (3.104),  (3.106),  and  (3.108)  for 
contaminant  transport  equations  and  in  the  following  equation  for  the  transport  equations  of  suspended 
sediments. 

(a)  when  ID  =  2 1 , 
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(a)  when  ID  =  22, 
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Subroutine  RINCS 

This  subroutine  is  to  read  in  chemistry  and  sediment  information. 

Subroutine  RINF_2 

This  subroutine  is  to  read  in  information  for  2-D  flow  simulations. 

Subroutine  RINIT_2 

This  subroutine  is  to  read  in  the  initial  condition  for  2-D  flow  and/or  transport  simulations. 

Subroutine  RINT_2 

This  subroutine  is  to  read  in  information  for  2-D  transport. 

Subroutine  SOL_BS 

This  subroutine  is  to  solve  the  transport  governing  equation  of  bed  sediments  at  the  NP-th 
river/stream  or  node. 

Subroutine  SOL_C2 

This  subroutine  is  to  solve  for  the  concentrations  of  dissolved  chemicals,  particulate  chemicals  on 
suspended  sediments,  and  particulate  sediments  on  bed  sediments  at  2-D  global  nodes.  This  subroutine 
serves  as  the  corrector  step  in  our  predictor-corrector  approach  to  solve  transport  equations  of  chemicals. 
The  Newton-Raphson  method  is  employed  to  solve  the  nonlinear  system  of  chemical  reactions.  The  full- 
pivoting  technique  is  incorporated  into  a  direct  solver  to  handle  linearized  matrix  equations. 

Subroutine  SOL_SS2 

This  subroutine  is  to  solve  the  concentrations  of  suspended  sediments  at  the  NP-th  node  in  the 
corrector  step,  where  the  time-implicit  scheme  is  employed. 

Subroutine  SOLVE 

This  subroutine  is  to  solve  band  matrix  equations  by  using  a  direct  solver.  The  band  matrix  equation 
is  triangulized  when  the  indicator  KKK  is  set  to  1,  while  back  substitution  is  implemented  when  KKK  =  2. 

Subroutine  STAR2M 

This  subroutine  is  to  determine  the  Lagrangian  values  (i.e.,  H,*,  u,*,  Vj*,  H2*,  u2*,  v2*,  H3*,  u3*,  v3* 
appeared  in  Section  3.2.1)  associated  with  each  characteristic. 

Subroutine  STORE_2 

This  subroutine  is  to  write  out  the  computer  results  at  desired  times  in  ascii  files  by  taking  the 
TECPLOT  format  for  achieving  the  purpose  of  post-processing  (e.g. ,  plotting).  Four  files  are  to  be  generated 
in  this  subroutine.  They  are  (1)  a  file  to  store  data  of  water  surface  elevation,  water  depth,  and  velocity,  (2) 
a  file  to  store  concentrations  in  water  body  and  interstitial  phase,  (2)  a  file  to  store  concentrations  of  both 
suspended  and  bed  sediments,  and  (4)  a  file  to  store  information  of  particulate  chemicals  adsorbed  on 
suspended  and  bed  sediments. 

Subroutine  STORE_2D 

This  subroutine  is  to  store  the  computer  results  at  desired  times  in  a  binary  file  for  further  plotting 
purposes. 


Subroutine  STRIP 
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This  subroutine  is  to  read  in  a  character  variable. 

Subroutine  SURE_2 

This  subroutine  is  to  analytically  determine  the  target  location  of  a  2-D  tracking  path  when  the 
single -velocity  approach  is  applied. 

Subroutine  SURE2M 

This  subroutine  is  to  analytically  determine  the  target  location  of  a  2-D  tracking  path  when  the 
single -velocity  approach  is  applied.  It  is  used  only  when  the  method  of  characteristics  is  employed  to  solve 
the  2-D  flow  equation. 

Subroutine  SURF_2 

This  subroutine  is  to  generate  boundary  arrays  for  implementing  boundary  conditions  in  2-D 
simulations. 

Subroutine  TBTABL 

This  subroutine  is  to  write  the  mass  balance  table  of  mobile  materials  at  desired  times. 

Subroutine  TRAK_21 

This  subroutine  is  to  determine  (1)  whether  or  not  a  particle  will  pass  through  the  working 
subelement  being  considered  and  (2)  the  location  of  the  target  location  if  the  particle  will  pass  through  the 
working  subelement.  This  subroutine  is  employed  when  the  starting  location  is  where  an  subelement  node 
is  located. 

Subroutine  TRAK_22 

This  subroutine  is  to  determine  (1)  whether  or  not  a  particle  will  pass  through  the  working 
subelement  being  considered  and  (2)  the  location  of  the  target  location  if  the  particle  will  pass  through  the 
working  subelement.  This  subroutine  is  employed  when  the  starting  location  is  not  where  an  subelement 
node  is  located. 

Subroutine  TRAK21M 

This  subroutine  is  to  determine  (1)  whether  or  not  a  particle  will  pass  through  the  working 
subelement  being  considered  and  (2)  the  location  of  the  target  location  if  the  particle  will  pass  through  the 
working  subelement.  This  subroutine  is  employed  when  the  starting  location  is  where  an  subelement  node 
is  located.  It  is  used  only  when  the  method  of  characteristics  is  employed  to  solve  the  2-D  flow  equation. 

Subroutine  TRAK22M 

This  subroutine  is  to  determine  (1)  whether  or  not  a  particle  will  pass  through  the  working 
subelement  being  considered  and  (2)  the  location  of  the  target  location  if  the  particle  will  pass  through  the 
working  subelement.  This  subroutine  is  employed  when  the  starting  location  is  not  where  an  subelement 
node  is  located.  It  is  used  only  when  the  method  of  characteristics  is  employed  to  solve  the  2-D  flow 
equation. 

Subroutine  TRANSP_2 

This  subroutine  serves  as  the  control  panel  to  compute  both  chemical  and  sediment  transport.  The 
predictor-corrector  numerical  approach  is  employed  to  solve  the  transport  equations  of  mobile  materials,  such 
as  dissolved  chemicals,  suspended  sediments,  and  particulate  chemicals  on  suspended  sediments.  During 
each  time  step,  The  predictor  step  is  first  implemented  to  obtain  the  intermediate  values  for  mobile  materials. 
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Sediment  equations  are  then  solved  by  using  the  Picard  to  handle  their  nonlinearity.  The  concentrations  of 
chemicals  are  the  last  to  solve  for.  The  set  of  nonlinear  ordinary  differential  equations  is  solved  with  the 
Newton-Raphson  method  node  by  node  in  the  corrector  step. 

Subroutine  TRILRN2M 

This  subroutine  is  to  generate  pointer  arrays  for  the  use  in  constructing  the  matrix  equation  for 
adjusting  2-D  flow  velocity  in  the  last  part  of  the  method  of  characteristic  approach. 

Subroutine  TSCONV 

This  subroutine  is  to  retrieve  desired  profiles  from  given  x-y  series  curves. 

Subroutine  TSCONVV 

This  subroutine  is  to  retrieve  desired  profiles  from  given  x-y  series  curves. 

Subroutine  TSCONV2 

This  subroutine  is  to  retrieve  desired  profiles  from  given  x-y  series  curves. 

Subroutine  TSFLOW 

This  subroutine  is  to  compute  the  data  needed  in  the  mass  balance  table  of  mobile  material. 

Subroutine  VALBDL_2 

This  subroutine  is  to  compute  the  velocity  at  a  specified  2-D  node  by  interpolation. 

Subroutine  VALBDL2M 

This  subroutine  is  to  compute  the  velocity  at  a  specified  2-D  node  by  interpolation.  It  is  used  only 
when  the  method  of  characteristics  is  employed  to  solve  the  2-D  flow  equation. 

Subroutine  WRKARY_2 

This  subroutine  is  to  prepare  working  arrays  for  particle  tracking  in  either  a  2-D  element  or  a  2-D 
subelement. 

Subroutine  WRKARY2M 

This  subroutine  is  to  prepare  working  arrays  for  particle  tracking  in  either  a  2-D  element  or  a  2-D 
subelement.  It  is  employed  only  when  the  method  of  characteristics  is  employed  to  solve  the  2-D  flow 
equation. 

Subroutine  XSI_2 

This  subroutine  is  to  determine  the  local  coordinates  of  a  specified  node  that  is  located  in  a  2-D 
quadrilateral  element,  given  the  original  Cartesian  coordinates  of  both  the  node  and  element  nodes. 
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